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Abstract 

The development of numerical methods for hyperbolic conservation laws has been a rapidly growing 
area for the last ten years. Many of the fundamental concepts and state-of-the-art developments can 
only be found in meeting proceedings or internal reports. This review paper attempts to give an 
overview and a unified formulation of a class of shock-capturing methods. Special emphasis will be on 
the construction of the basic nonlinear scalar second-order schemes and the methods of extending these 
nonlinear scalar schemes to nonlinear systems via the exact Riemann solver, approximate Riemann 
solvers, and flux-vector splitting approaches. Generalization of these methods to efficiently include 
real gases and large systems of nonequilibrium flows will be discussed. The performance of some of 
these schemes is illustrated by numerical examples for one-, two- and three-dimensional gas-dynamics 
problems. 
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I. Introduction 


This paper attempts to give an overview and a unified formulation of a class of shock-capturing 
methods. Much of the mathematical theory is omitted. However, when appropriate, sufficient refer- 
ences will be provided. It is assumed that the reader is familiar with the fundamentals of numerical 
analysis. Before getting into the main discussion, some pertinent terminology will be covered. The 
specifics to be addressed, and some aspects of numerical schemes for hyperbolic conservation laws will 
be summarized 

Terminology : In this paper, the terms explicit and implicit schemes refer to time discretization, whereas 
the terms symmetric and upwind schemes refer to spatial discretization. The order of accuracy for 
time-accurate calculations refers to both the time and spatial discretization. On the other hand, the 
order of accuracy for steady-state calculations most often refers to the spatial discretization only. 

Upwind-differencing schemes attempt to discretize hyperbolic partial differential equations (PDE’s) 
by using differences biased in the direction determined by the sign of the characteristic speed. Symmet- 
ric or centered schemes, on the other hand, try to discretize hyperbolic PDE’s without any knowledge 
of the sign of the characteristic speed. 

Shock-capturing schemes tend to treat shocks as a continuum, as opposed to shock-fitting, where 
almost always a priori knowledge of the shock location is needed. 

Specifics to be Addressed: In this paper, only conservative finite-difference methods for hyperbolic 
conservation laws (i.e., for inviscid compressible flows) are addressed. The formulations are Eulerian, 
and the main emphasis is state-of-the-art of a class of high-resolution shock-capturing methods for the 
last ten years. Only initial-value problems are considered. For compressible Navier-Stokes calculations, 
the physical problems considered here are assumed to be inviscid-dominated in the sense that moderate 
or strong viscous shock waves are present in the flow field such that high-resolution shock-capturing 
techniques are required. Thus, the numerical procedure proposed here for Navier-Stokes calculations 
is that a second-order, central-difference approximation is used for the diffusion terms and a high- 
resolution shock-capturing method is used for the inviscid part of the Navier-Stokes equations (Euler 
equations). 

Hierarchy of Conservative Schemes for Hyperbolic Co nservation Laws : The Hierarchy of conservative 
schemes for hyperbolic conservation laws can best be illustrated by figure 1.1. Let St be the set of 
all existing conservative schemes of any order for hyperbolic conservation laws which is the entire 
region shown in figure 1.1. We can break St into two parts, Sup and Sc, where Sup is the set 
of all existing upwind schemes of any order. Furthermore, let Seno be the set of all essentially 
nonoscillatory (ENO) [1-3] schemes of any order, let Stvd be the set of all total variation diminishing 
(TVD) schemes [4-12] of any order and let S M be the set of all monotone schemes [13-14]. Then 
Sm C Stvd c Seno c St- In other words, the set of monotone schemes is the smallest set and is a 
subset of the set of TVD schemes. The set of all TVD schemes in turn is a subset of the ENO schemes. 
Definition and properties of these schemes will be described inside the text. This paper covers only a 
small subset of the shock-capturing schemes, namely, the TVD schemes. 

Cl assical vs. Modern Shock- Capturing Methods : From an historical point of view, shock-capturing 
methods can be classified into two general categories: namely, classical and modern shock-capturing 
methods. In the case of classical shock-capturing methods, numerical dissipation terms are either linear 
such that the same amount is applied at all grid points or contain of adjustable parameters [15-17]. 
Classical shock-capturing methods only exhibit accurate results for smooth or weak shock solutions 
and are not robust enough for strong shock wave calculations. For strong shock waves, classical shock- 
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capturing methods either result in oscillatory solutions across the discontinuities and/or nonlinear 
instability. For modern shock-capturing methods, however, the numerical dissipations are nonlinear 
such that the amount varies from one grid point to another and usually consist of automatic feedback 
mechanisms to control the amount of numerical dissipation. These schemes are stable for nonlinear 
scalar hyperbolic conservation laws [1-14] . Numerical dissipation terms similar to those of Jameson et 
al. [18] seem to fall in between the classical and modern shock-capturing methods. 

Applicability of Scalar Schemes to Systems of Hyperbolic Conservation L aws : Basic modern shock- 

capturing methods for hyperbolic conservation laws are developed for linear and/or nonlinear scalar 
hyperbolic conservation laws. Extension of scalar methods to nonlinear systems is accomplished by as- 
suming certain physical models or by local linearization. The mathematical foundation relies mainly on 
the scalar case. There is no identical theory for nonlinear systems or the multidimensional counterpart. 
These schemes are formally extended to one- or higher-dimensional nonlinear systems of hyperbolic 
conservation laws via the so called Riemann solvers (to be defined in section 4) and are evaluated 
by numerical experiments. Based on these facts, a major part of the discussion will be on modern 
shock-capturing schemes for the nonlinear scalar case and the methods of extending these nonlinear 
scalar schemes to nonlinear systems via the different types of Riemann solvers. However, numerical 
examples for one- and higher-dimensional nonlinear systems in gas dynamics will be stressed. 

The discussion on the developments of conservative difference schemes for hyperbolic conservation 
laws is the author’s personal interpretation. The illustrations of the types of schemes and the areas 
of applications reflect the author’s experiences and preferences for certain schemes. No attempt has 
been made to present a unified comparison except for the one-dimensional shock-tube problem. 

Outline of Paper The outline of the paper is shown in the table of contents. First, the basic properties 
of hyperbolic conservation laws and several schemes for linear hyperbolic equations will be reviewed. 
Then the various aspects of shock-capturing schemes for nonlinear scalar hyperbolic conservation laws 
will be discussed. These include monotone and first-order upwind schemes, deficiency of classical 
shock-capturing schemes, and methods of extending first-order TVD schemes to higher-order. Since 
the general theory of ENO schemes is very involved and still in the development stage, it will not be 
discussed here. In section IV, formal extensions of nonlinear scalar TVD schemes to one-dimensional 
nonlinear systems of hyperbolic conservation laws will be reviewed. A method of extension, which is 
widely used and practical in terms of complex fluid dynamics problems, will be stressed. Time-accurate 
as well as steady-state calculations for one-, two- and three-dimensional practical applications will be 
illustrated when appropriate. 


5 



II. Preliminaries 


The main difficulty in solving hyperbolic PDE’s is the need to allow for discontinuous solutions 
even when the initial data are smooth. For constant-coefficient hyperbolic PDE’s, well known stable, 
finite-difference methods are available in standard textbooks; see for example: Richtmyer and Morton 
[19], Mitchell [20], and Garabedian [21]. The theory is more complex for nonlinear hyperbolic PDE’s. 
In order to motivate the ideas and set up the basic notations for nonlinear hyperbolic PDE’s, some of 
the schemes originally designed for constant-coefficient cases are reviewed in this section. 

2.1. An Upwind Scheme for Linear Hyperbolic PDE’s 
Consider the constant-coefficient scalar hyperbolic PDE 


du du 
di + a d~x 


( 2 . 1 ) 


where a is a real constant. Let u ” be the numerical approximation of the solution of (2.1) at x = j Ax 
and t = nAt , with Ax the spatial mesh size and At the time-step. Let A = According to the 
characteristic theory of hyperbolic PDE’s (direction of wave propagation), a simple first-order accurate 
explicit upwind scheme can be written as 


u 


n+l 


Introducing the notation 


a+ = - (a + |a|); 


- u n 
3 

a < 0 

u n 
U 3~ 1 

a > 0 

1 , 

- M), 

= 2 ( “ 


the scheme for positive or negative a can be rewritten as one equation 


( 2 . 2 ) 

(2.3) 


n+l n 


= Uj - A [o + (u? - «“_,) + «-(«”+ 1 - «*?)]. 


(2.4) 


Using the relationship between a^, a and |a|, the scheme again can be written as 


n+l 


= «? - .»(< 


3+1 




I )+-HK +1 -2«” + uj.i). 


(2.5) 


Most often one recognizes the first form (2.4) as an upwind scheme but the second form (2.5) is less 
obvious. In this paper, the second form is preferred because when one extends the scheme to nonlinear 
equations and systems of equations, the second form is more compact to discuss and more efficient 
in terms of operations count [22-24]. Higher-order upwind schemes can be obtained by increasing the 
stencils of the first-order scheme in the appropriate upwind directions. 


Consider a system of hyperbolic equations 


dU dU n 
+ - 0, 


( 2 . 6 ) 


where U is a vector with m elements and A is an m x m constant matrix with real eigenvalues. Let 
W — R~ l U and R~ 1 AR = A. One can transform the above system to a diagonal form 


-T-+A— =0, A - diag(a 2 ), / = 1, ...,m. (2.7) 

ot ox 

Here diag(a*) denotes a diagonal matrix with diagonal elements a 1 . Applying the scalar scheme to the 
new variables, one obtains a scheme for the system case: 
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(2.8a) 


W? +1 = w i - - W"_ t ) + ||A|(W" +1 - 2 W? + W^), 

with 

|A| = diag(|a*|). (2.8b) 

This form looks exactly like the scalar case except it consists of m equations. Transforming back to 
the original variables, the scheme takes the form 


y” +1 = t? - |a(h; + 1 - £/;_,) + ||A|(ti; +1 - 21/; + u^). 

( 2 . 9 a) 

with 

|A| - ii|A|i 2 -1 , 

( 2 . 9 b) 

or 


U " +1 = ti; - a[a+(h; - £/;„,) + a-(v; +1 - 17)], 

(2.10a) 

where 


A+= 1 (A+|A|), 

(2.10b) 

a-=\(a-\a\). 

(2.10c) 

Introducing a new notation “F J + 1”, consider a scheme of the form 



(2.11) 

where P,- . 1 is sometimes referred to as a “numerical flux function” . Note that 

3 "r 2 

will be used throughout this paper. For the previous example ( 2 . 10 ) 

the notation Fj + i 

F jJth = \[A(V, +i + Vj) - \A\(U, +i - Vj)]. 

(2.12) 


2.2. Centered (Symmetric) Schemes for Hyperbolic PDE’s 

Several popular, spatially centered, second-order accurate schemes for the constant-coefficient scalar 
hyperbolic PDE’s are as follows: 

(i) Crank-Nicholson method: 


3 1 2 v •?+ 

(ii) Lax-Wendroff method: 

Aa 


n+i + ^( u n+l _ u n+U = n _ — 

3 1 O v j+1 J- J 9 ' 3 + 1 3~ 1/* 


1 = «? - yK+i - «J-i) + 5 AV(«j + 1 - 2« j + 


( 2 . 13 ) 


( 2 . 14 ) 


(iii) MacCormack method: 


u 3 1} = " Aa («"+i “ u ? )’ 


( 2 . 15 a) 
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(2.15b) 



Note that the Lax-Wendroff method can be rewritten as 

u? +1 = u"-A (2.16a) 

with the numerical flux function h-, i 

J+ 2 

h J + i = - [a(u i+1 + Uj) - X a 2 (u j+1 - tiy)]. (2.16b) 

The Crank-Nicholson method can be expressed similarly. 
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III. Schemes for Nonlinear Scalar Hyperbolic Conservation Laws 

The main theory for modern shock-capturing methods for nonlinear hyperbolic conservation laws 
considered in this paper relies heavily on the basic first-order upwind and the Lax-Wendroff methods. 
However, the resulting higher-order modern shock-capturing methods, which are designed for the 
nonlinear hyperbolic conservation laws, when applied to constant-coefficient cases, are no longer linear 
finite-difference schemes (i.e. , they are truly nonlinear finite-difference methods). This fact will be 
stressed in the current section. 

3.1 Conservative Schemes and a Shock-Capturing Theory 

The stability analysis of difference schemes for linear hyperbolic PDE’s is very well established. 
However, the stability analysis of difference schemes for nonlinear hyperbolic PDE’s is less developed. 
In general, the stability theory for linear difference equations is of use in checking the “local” stability of 
linearized equations obtained from truly nonlinear equations. However, in many instances when strong 
discontinuities are present, local stability is neither necessary nor sufficient for the nonlinear problems. 
One traditional remedy for removing instabilities is to introduce a “linear” numerical dissipation (or 
“artificial viscosity” or “smoothing term”) into the difference schemes. One can do so by designing 
the scheme to be dissipative [19]. 

The majority of practical applications in fluid dynamics for the last decade is based on this traditional 
approach of adding an additional dissipation term to the numerical scheme to improve nonlinear 
instability. However, this approach alone will not guarantee convergence to a physically correct solution 
in the nonlinear case. 

Lax and Wendroff [25] showed that the limit solution of any finite-difference scheme in a conservation 
form which is consistent with the conservation laws satisfies the jump conditions across a disconti- 
nuity automatically. This was a conceptual breakthrough which enabled the direct discretization of 
the conservation laws by introducing the notion of numerical dissipation. However, weak solutions 
(solutions with shocks and contact discontinuities) of hyperbolic conservation laws are not uniquely 
determined by their initial values; an entropy condition is needed to pick out the physically relevant 
solutions. The question arises whether finite-difference approximations converge to this particular 
solution. It is shown in references [13,14] that in the case of a single conservation law, monotone 
schemes (to be defined later) always converge to the physically relevant solution. If the scheme is not 
monotone, then it must be consistent with an entropy inequality for the assurance of convergence to 
a physically relevant solution [26,27]. Thus monotone schemes possess many desirable properties for 
the calculation of discontinuous solutions. Moreover, first-order upwind schemes share most of the 
properties of monotone schemes. The following is an introduction to monotone and first-order upwind 
schemes. 


3.2. Monotone and First-Order Explicit TVD Schemes 
Consider the scalar hyperbolic conservation law 

du df(u) 


dt + dx 


0 , 


(3.1) 


where a(u) = df jdu is the characteristic speed. A general three-point explicit difference scheme in 
conservation form can be written as 


„ ? +> = u »-a(a“ + ± 


(3.2) 
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where h” +i = , u” +1 ). The numerical flux function is required to be consistent with the 

conservation law in the following sense 


Rewrite equation (3.2) as 


h(u 3 ,uj) = f(u 3 ). 


n-f-1 

3 




(3.3) 


(3.4) 


The numerical scheme (3.4) is said to be monotone if G is a monotonic increasing function of each of 
its arguments. Monotone schemes produce smooth transitions near discontinuities, but they are only 
first-order accurate [13,14]. Examples of monotone schemes are the Lax-Friedricks scheme [14], the 
Godunov scheme [28], and the Engquist and Osher scheme [29]. 

In general, the class of first-order upwind schemes is larger than the class of monotone schemes. 
Not all first-order upwind schemes are monotone schemes. For example, the Godunov method is a 
first-order monotone upwind scheme and is of the form 



min Uy <u<u J+ i/( u ) u j ^ «y+i 

max Uy > u > Uy+1 /(u) u 3 > 


(3.5) 


The Engquist and Osher method is again a first-order monotone upwind scheme. However, the first- 
order upwind schemes of Huang [30] and Roe [31] are not monotone schemes. All of these popular 
first-order upwind schemes (with explicit Euler time-discretization) can be cast in the following form: 


w ” +1 = = - A D? - \D%. 

Here, D i represents some forward difference of the / or u. For example, D\ can be 

Di — /j+i) Uj, Uj+i)(/j+i ~ fj) 


(3.6) 


(3.7a) 


or 

Di — A 2 {fj) />+i5 u ji u .?+i)( w y+i — (3.7b) 

and D 2 represents some backward difference of the / or u. For example, D 2 can be 

D 2 — Si(/j-i) fj, u j—i , — fj- 1 ) (3.8a) 


or 

D 2 = B>2 {fj- 1, fj, Uj- 1 , Uj)(uj - tty _ 1 ) . (3.8b) 

Here A 1 , A 2 , B 1 , and B 2 are some known functions of the arguments indicated above. As an example, 
consider the Engquist and Osher scheme, where the D\ and D 2 for any convex flux function / are 


= D, = (3.9a) 

with 

/+ = /(max(uj, u)), 
fj = /(min(«j,u)), 
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(3.9b) 

(3.9c) 



and u is the sonic point of /(ti); i.e., /'(u) = 0. 

In [30], Huang introduced a first-order accurate upwind scheme 

u " +1 = «" - 5 [» - •8»(«“ +1 )](/?+. - /”) - * [i + •*»(«?_*)](/; - /;_i). ( 3 -io) 

She was vague in defining &j+i for a general flux function /, but for Burgers’ equation, she explicitly 
defined 


Here 


a ) i i - (°j'+i + 


°i = » t 1 - s 8“(°)+i)](/j+i “ fi)> 


which has the same form as (3.7a), and 


02 = i[l + sgn(o 3 _i)](/,' - /j-i), 


(3.11) 

(3.12) 

(3.13) 


which has the same form as (3.8a). 
In [31], Roe defined 


°i+± - 


_ f (/j+i - /j)/ a j+ A j+|“ ^ 0 


A j+i“ = °- 


\ <*K) 

where i« = + i — u^ . This is equivalent to Huang’s method for Burgers’ equation. 


(3.14) 


With the definition (3.14), scheme (3.10) can be rewritten as a three-point central difference method 
plus a numerical viscosity term 


."+1 — ..n 


«“ ~ olf" + 1 - - K + il A i + i«" + l«?-il A ,--i« n ] 


(3.15) 


Now, using definition (3.14), D i = | [ a j+% ~ l a j+i|]( u i+i — *fj), which has the same form as (3.7b), 
and D 2 — \ [°j-i + \ a j- 1 1] { u j ~ ttj-i), which has the same form as (3.8b). 


The numerical flux function as a function of D\ can be written as 

h 3+\ = fj + 

Or, one can express (3.15) as (3.2) with 

h i+i = \[f> + fi+i - lK<*,-+i) A ,-+i“]> 


(3.16) 


(3.17a) 


and 


V'K+i) = K+i|- (3.17b) 

tj) is sometimes known as the coefficient of the numerical viscosity term. In this paper, I prefer to 
use equation (3.17a) as the form of the first-order upwind numerical flux function. This form of the 
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numerical flux function (3.17a) is not a common notation. But we can see later that representation 
(3.17a) is quite useful for the development of second-order TVD schemes, especially via the modified- 
flux approach [6]. Moreover, (3.17a) provides a more compact form for extension to nonlinear systems 
[22]. 


It is well known that (3.10) or (3.15) is not consistent with an entropy inequality, and the scheme 
might converge to a nonphysical solution. A slight modification of the coefficient of numerical viscosity 
term [6], 


v-(*) 


\z\ \z\> Si 

(z 2 + Si)/2Si \z\ < Si, 


(3.18) 


can remedy the entropy violating problem. Here ip(z) is an entropy correction to |z|, where Si is a 
small positive parameter or a function of z (see reference [32] for a formula for 6 1 ). Other ways of 
modifying (3.17b) to satisfy an entropy inequality can be found in [33-35]. One can view the size of 
as a measure of the amount of numerical dissipation for the first-order upwind numerical flux (3.17a). 

= 0 is the least dissipative, and the larger the the more dissipative the scheme becomes. Section 
5.7 discusses the use of for steady-state, blunt-body hypersonic flows. 

If we define 


C ± (z) = \[^z)±z ] |, 


then, this upwind scheme can be written as 


u 


n+ 1 


= «? + XC~(a^)A j+ ^u n - 


(3.19) 


(3.20) 


In other words, this conservative scheme can be viewed as a generalization of the Courant et al. 
nonconservative upwind scheme [36]. 

3.3. Deficiency of Classical Shock-Capturing Methods 

Although monotone schemes possess many desirable properties for the calculation of discontinuous 
solutions, they are only first-order accurate [13,14]. For complex flow-field structures, monotone and 
first-order upwind schemes are too diffusive. They cannot produce accurate solutions for complicated 
flow fields with a reasonable grid spacing. One needs higher-order shock-capturing methods. In the 
last ten years the emphasis has been on the development of better methods for problems with shocks. 
As discussed in the introduction, one can loosely divide higher-order shock-capturing methods into two 
classes. The classical one uses linear numerical dissipation; i.e., it uses the same amount everywhere 
or consists of adjustable parameters for each problem. The modern one uses nonlinear numerical 
dissipation. The amount varies from grid point to grid point and is built into the scheme with hardly 
any adjustable parameters. 

Higher-order accurate classical upwind and symmetric shock-capturing schemes suffer from the 
following deficiencies: (l) they produce spurious oscillation whenever the solution is not smooth, (2) 
they may develop nonlinear instabilities when discontinuities are encountered, and (3) the scheme 
may select a nonphysical solution. Figure (3.1) shows an example of spurious oscillations associated 
with the classical shock-capturing method where the Burgers’ equation is solved by the Lax-Wendroff 
method. Here the flux function /(u) = u 2 / 2. The initial condition is taken to be a sine wave and the 
boundary condition is taken to be periodic. The solid lines are the exact solutions at two different 
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times and the circles are the solutions computed by the Lax-Wendroff method. At the time when the 
solution is still smooth the computed solution matches with the exact solution very well. However, 
at the time when the solution has developed into a shock, the scheme produces oscillations across the 
shock. The oscillations near the shock remain the same as the mesh is refined. 

There are two classes of modern shock-capturing schemes which are more appropriate for the cal- 
culation of weak solutions, namely the TVD and ENO schemes. In addition, these schemes should 
be consistent with an entropy inequality, and second- or higher-order in smooth regions, and should 
produce high-resolution solutions near shock and contact discontinuities. Most of the available higher- 
order TVD and ENO schemes possess these properties. The main distinction between ENO and TVD 
methods is that ENO schemes can retain the same spatial order of accuracy even at points of extrema, 
whereas TVD schemes reduce to spatially first-order at these locations. TVD schemes are a subset 
of ENO schemes and are more efficient in terms of operations count. ENO schemes are still in the 
development stage, whereas TVD schemes are more established in the sense of application to a wide 
range of multidimensional gas-dynamics problems. Only TVD schemes will be discussed here. Before 
a detailed discussion of the theory and the method of construction, the performance of a second-order 
TVD schemes will be illustrated. A second-order TVD scheme developed by Harten was* applied to 
the same Burgers’ equation at the same two time instances as the Lax-Wendroff schemes as shown in 
figure (3.2). The solutions are very smooth near the shock. 

The next section is devoted to the introduction of TVD schemes. The definition and sufficient 
conditions for a scheme to be TVD will first be covered and then some first-order TVD scheme examples 
will be given. It turns out that all the monotone and first-order upwind schemes are first-order 
TVD schemes, and all first-order TVD schemes generate monotonic shock profiles. Unlike monotone 
schemes, not all TVD schemes are automatically consistent with an entropy inequality. Consequently, 
some mechanism may have to be explicitly added to TVD schemes to enforce the selection of the 
physical solution. An example is the entropy correction t p(z) ((3.18)) to \z\ for the first-order Roe 
scheme. It is emphasized here that the TVD property is only valid for homogeneous scalar hyperbolic 
conservation laws. For nonhomogeneous hyperbolic conservation laws, in order for the source term 
to not influence the TVD property, it is restricted to a special class of functions. For example if the 
source term is contractive in the sense of stiff ordinary differential equations, it is expected that the 
source term will not influence the TVD property. 

3.4. TVD Schemes: Background 

Consider a one-parameter family of five-point difference schemes in conservation form, 


n+ 1 


+ X6(h 


n+ 1 

J+ 2 


J ' 2 


A(1 -»)(*?+* 


(3.21) 


where 0 < 6 < 1, h" L = and = h(< 


n+1 n+1 n+ 1 n+1 


<+*). The same 


n+l 

"j'+£ ~3 ’ ~J + 1’ -J + Z'’ j+A , *t u 'j-l’“j ’ “j + 1 ’ “j + 2 

numerical flux function with a different time-index appears on both sides of the equation. Let 


h. 


n + l 
0+ 2 


be another numerical flux function. Then (3.21) can be rewritten as 


(3.22) 


u n+1 = u n 
3 J 


A(/i J + i - h 7 _i). 


(3.23) 


Here = h(uy_ 1 ,u”,u" +1 ,u” +2 ,u"^ 1 1 ,u” +1 ,u”^ 1 1 ,u”^ 2 1 ) and is consistent with the conservation 

law (3.1) in the following sense 
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/i(tt,u,u,u,u,u,ti,u) = /(u). (3.24) 

This one-parameter family of schemes contains implicit as well as explicit schemes. When 9 = 0, 
(3.21) is an explicit method. When 9 ^ 0, (3.21) is an implicit scheme. For example, if 9 = 1/2, the 
time-differencing is the trapezoidal formula, and if 9 — 1, the time-differencing is the backward Euler 
method. Other forms of difference schemes can also be analyzed. However, for implicit methods it 
will be more difficult to analyze the TVD property. To simplify the notation, rewrite equation (3.21) 
as 


L • u n+1 = R - u n , 

where L and R are the following finite-difference operators: 


(3.25) 


( L ■ u )y = u 3 + - h 3-0 

(R • u) 3 = Uj - A(1 - B)(h j+ 1 - h,-_ i). 


The total variation of a mesh function u n is defined to be 

CO oo 

TV(u”)= J2 K +1 - U "l= £ l A i+i“"l 


(3.26a) 

(3.26b) 


(3.27) 


J = — oo 


3 =- oo 


Here the general notation convention 


A J+i Z=Z 3 + l -Z > ( 3 - 28 ) 

for any mesh function z is used. The numerical scheme (3.21) for an initial-value problem of (3.1) is 
said to be TVD if 


TV (u n+1 ) < TV(u n ). (3.29) 

The following sufficient conditions for (3.21) to be a TVD scheme are due to Harten [6]: 

TV (R - u n ) < TV (u n ) (3.30a) 

and 

TV(L • u n+1 ) > TV{u n+1 ). (3.30b) 

Assuming that the numerical flux h J+ 1 in (3.21) is Lipschitz continuous, (3.21) can be written as 


n+ 1 


-X9 


C . i A • , iu 

3+± J+2 


c 


n-f 1 


i A ,■-*«) =«j+a(i-^) 


C y+i A j+$ u 


i A_iu 
3 ~ 2 3 2 


(3.31) 


where CT , = <5^ (u, T i, Uj-tj, u, ±2 ) are some bounded functions. Then Harten further showed 

3^- 2 

that sufficient conditions for (3.30) are 
(a) if for all j 
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(3.32a) 


> 0 

+ C- + j<l, (3.32b) 


< 0 (3.33) 

for some finite C. For example, when 0 = 0 and = C ± = ± z\ as defined in (3.19), the 

resulting scheme (3.21) is a first-order explicit upwind scheme, whereas when 9 = 1 with the same C^, 
the scheme is a first-order implicit upwind method. Both of these methods satisfy conditions (3.32) 
and (3.33). By examining all the first-order upwind schemes in section 3.3, it can be shown that they 
are first-order TVD schemes. Conditions (3.32) and (3.33) are very useful in guiding the construction 
of higher-order-accurate TVD schemes which do not exhibit the spurious oscillations associated with 
the more classical second-order schemes. Other necessary and/or sufficient conditions for semi-discrete 
difference methods for nonlinear hyperbolic PDE’s can be found in references [37,12]. 


and 

(b) if for all j 


C t +i = A (1 - »)cf +i 


C+ +i +C^=X(l-t)lC+ +i 


-oo < C < -AtfCf, , 


3.5. Higher-Order TVD Schemes 

The author is aware of primarily four different (and yet not totally distinct) design principles for the 
construction of high-resolution TVD schemes. They are (1) hybrid schemes such as the flux-corrected 
transport (FCT) of Boris and Book [38], Harten [39], and van Leer [40]; (2) second-order extension 
of Godunov’s scheme by van Leer [4], and Colella and Woodward [5]; (3) the modified-flux approach 
of Harten [6]; and (4) the numerical fluctuation approach of Roe [7] and Sweby [41]. Also, Osher [42] 
has subsequently extended the first-order scheme of Engquist and Osher to second-order accuracy by 
using the above ideas. More recently, Jameson and Lax [12] extended the TVD idea for multi-point 
schemes. The following is a subjective interpretation of these design principles. 

(1) The flux-corrected transport scheme is a two-step hybrid scheme consisting of a combined first- 
and second-order scheme. It computes a provisional update from a first-order scheme, and then filters 
the second-order corrections by the use of flux limiters to prevent occurrence of new extrema. 

The idea of the hybrid scheme of Harten or van Leer is to take a high-order- accurate scheme 
and to switch it explicitly into a monotone first-order-accurate scheme when extreme points and 
discontinuities are encountered. 

(2) van Leer observed that one can obtain second-order accuracy in Godunov’s scheme by replacing 
the piecewise-constant initial data of the Riemann problem with piecewise-linear initial data. The slope 
of the piecewise-linear initial data is chosen so that no spurious oscillations can occur. Woodward and 
Colella further refined van Leer’s idea by using piecewise-parabolic initial data. 

(3) The modified-flux TVD scheme is a technique to design a second-order accurate TVD scheme 
by starting with a first-order TVD scheme and applying it to a modified flux. The modified flux is 
chosen so that the scheme is second-order in regions of smoothness and first-order at points of extrema. 
Details of the construction of this scheme can be found in reference [6]. A discussion will be presented 
in a later section. 
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(4) The numerical fluctuation approach of Roe is a variation of the Lax-Wendroff scheme. Roe’s 
variation depends on an average function. The average function is constructed (in such a way that 
spurious oscillations will not occur) by the use of flux limiters. As a matter of fact, under certain 
assumptions, a form of Roe’s second-order scheme is equivalent to the modified-flux approach. 

Most of the above methods can also be viewed as three-point central-difference schemes with an 
“appropriate” numerical dissipation or smoothing mechanism. “Appropriate” here means automatic 
feedback mechanism to control the amount of numerical dissipation, unlike the numerical dissipation 
used in linear theory. Design principles (2)-(4) are more closely related to the mathematical notion of 
TVD schemes devised by Harten. Hybrid types of schemes similar to design principle (1) do not fit in 
the same mathematical notion and will not be discussed in this paper. 

In general, the basic idea of the above design principles is to construct a higher-order scheme which 
has properties similar to a first-order TVD scheme so that spurious oscillations cannot be generated. 
The main mechanisms for satisfying higher-order TVD sufficient conditions are some kind of limiting 
procedures called limiters (or flux limiters). They impose constraints on the gradients of the dependent 
variable u or the flux function /. For constant coefficients, the two types of limiters are equivalent. 
One can obtain a second-order TVD scheme by modifying an upwind scheme or a centered scheme 
with proper limiters; i.e., if the scheme so constructed satisfies the TVD sufficient conditions. For 
nomenclature purposes, the term “upwind” or “symmetric” TVD schemes will be used to denote the 
original scheme before the application of limiters. Another way of distinguishing an upwind from 
a symmetric TVD scheme is that the numerical dissipation term corresponding to an upwind TVD 
scheme is upwind- weighted, whereas the numerical dissipation term corresponding to a symmetric 
TVD scheme is centered. This generic use of the notion upwind and symmetric TVD schemes no 
longer has its traditional upwinding and centering meanings. In general, symmetric TVD schemes are 
simpler than the upwind TVD schemes. This point will become apparent later. 


3.5.1. Higher- Order Upwind TVD Scheme 

Now we turn to discuss the derivation of higher-order schemes in conjunction with TVD properties 
(i.e., the method of obtaining higher-order schemes using limiters). There are many variations but 
basically they can be divided into two general approaches. The so-called “MUSCL” (monotonic 
upstream schemes for conservation laws) approach due to van Leer and the others, hereafter grouped 
under the non-MUSCL approach. The methods of van Leer [4], Harten [6], Roe [7], Sweby [41], and 
Osher and Chakravarthy [11] for upwind TVD schemes, and the methods of Davis [8], Roe [9] and Yee 
[10] for symmetric TVD schemes are typical examples of the MUSCL and non-MUSCL approaches. 
Note that one can obtain higher-order upwind TVD schemes via either the MUSCL or non-MUSCL 
approach, whereas one can obtain symmetric TVD schemes only via the non-MUSCL approach. 

For simplicity, take the forward Euler time-differencing of the form 


n-j- 1 


3 + i 



(3.34a) 


with 

Vff = g [^' +1 + + ^3+^]' (3.34b) 

Here the higher-order numerical flux ± is used to distinguish it from a first-order numerical flux 

h-, i. Also the numerical fluxes described below will be a first-order time-discretization for the 

3 ~ r 2 

MUSCL and Osher and Chakravarthy schemes. One would not recommend the use of the explicit 
Euler time-discretization for these two methods, since if the limiters are not present, linear stability 
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analysis shows that these two methods are unconditionally unstable. Forms that are second-order in 
time will be discussed at the appropriate places. 

MUSCL (Monotonic Upstream Schemes for Conservation Laws) Approach : Consider a three-point 

explicit difference scheme (2.2) in conservation form, 


= «? - A[k(«?,«J +l ) - (3 35) 

where the numerical flux h J + x is a function of Uj and tey+i. Use the short-hand notation A J+ 1 = 
Uj+i — uy (i.e., deleting the u from A J+ ±u; A J + i and A will be used interchangebly throughout 
the text) and consider a first-order upwind numerical flux function of Roe, 

h j+± ~ 2 [h + /i+i “ K+£I A j+£ • (3.36) 

In reference [4], van Leer observed that one can obtain spatially second-order accuracy in Godunov’s 
scheme by replacing the piecewise constant initial data of the Riemann problem with piecewise linear 
initial data. The MUSCL approach applied to the above first-order numerical flux (instead of Go- 
dunov’s scheme) to obtain a spatially higher-order differencing is to replace the arguments Uj+ \ and 
Uj by u R k and u L where u R and u L are defined as follow: 

3t j 3 i 2 

U f +1 = U 3+1 - \ [(! - »?) A J+ f + (! + > (3.37a) 

u y + 1 = u 3 + 4 [(! “ T ?) A i-f + (1 + *7)Ay + i] • (3.37b) 

Here the spatial order of accuracy is determined by the value of r?: 

rf = — 1, fully upwind scheme 
rj = 0, Fromm scheme 

r t = 1/3, third-order upwind-biased scheme 
rj — 1, three-point central-difference scheme 

Various “slope” limiters are used to eliminate unwanted oscillations. A popular one is the “minmod” 
limiter; it modifies the upwind-biased interpolation as follows: 

u f+i = u 3+1 - i[(l - rj) A j+ i + (l + r^A^i], (3.38a) 

= Uj + -[(1 - y)Aj _ i + (1 + y)Aj + i], (3.38b) 


where 


A i+i = minmod (Ay + i,wA 3 _i), 
A ,+i = minmod (A i+ 1 ,a>A 3+ |), 


minmod(a:,a;y) = sgn(z) - max] 0,min |z|,a>ysgn(:r) 


(3.38c) 

(3.38d) 

(3.38e) 
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and 1 < uj < with r) ^ 1. Therefore, a spatially higher-order scheme can be obtained by simply 
redefining the arguments of h J± i ; i.e. 


,n+i — 


tiJ-A 


h(uf + 1 , u L i+ 1 ) - /i(uf_ i , i ) 


j+ 


(3.39) 


Applying the above to the first-order Roe scheme, the second-order numerical flux by the MUSCL 

approach denoted by h VR k is 

3' 2 


*ri*=*i «?**.«* 


3 + 


. R 


3+ i 5 3+ ■ 


/(«f+») + /(«?+*)- 


J+: 


/(•&*)-/(«?+*) 




(*?+*- «f + t) 


J+- 


For n = — X, A„_l i = A„__i can be the limiter function as follows: 

* 1 J T 2 J 2 


A ;+§ “ | A j-|[( A y+i) 2 + *»] + A, + *[(A J . i ) 2 + «|/ 


(3.40) 


( A i+i) 2 + ( A i-i) 2 + 2 ^2 


(3.41) 


The parameter 10 7 < S 2 < 10 5 is a commonly used value in practical calculations. 

One way to obtain a second-order time-discretization (in addition to higher-order spatial discretiza- 
tion ) is to replace the forward Euler time-discretization by some linear multistep method or by the 
Runge-Kutta-type of time-discretization. Another way is to redefine u R i and u L k by the following: 

3 r 2 3 i 2 


< + i =(«" + * + \h) ( 3 - 42a ) 

«?+* =(«£,* 5 - lfe+i). (3.42b) 

Here gj is the same as A J+ i in equations (3.41) or the minmod function (3.38e) (with u; = 1, x 

n-)- 1 

replaced by A J + i and y replaced by Aj_i). The quantity u } 2 is 

U T 1 = u " - \ { f [(“? + |?i)]-/ [(«? - } • ( 3 - 43 ) 

Modified- Flux Approach ( Harten ): Now, the second-order TVD scheme of Harten [6] is considered. 
His method is sometimes referred to as the modified-flux approach. Apply the first-order TVD scheme 
to an appropriately modified flux 


/«—(/ + »)• (3.44) 

The second-order numerical flux looks exactly like the first-order scheme, except it is a function of 
/ = (/ 4- g) instead of /. Thus, the second-order numerical flux denoted by h H L is 

3* o 


ftf+i = \[f, + fi+i - ^( s ,-+i) A )+iL 


with 



(3.45a) 


(3.45b) 
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Vi =°)+i+ij+i' 


(3.45c) 


(ffj+i-ffy)/Vi A )+i^° 

0 i, + f=0, 


(3.45d) 


where the function i p(z) is an entropy correction to \z\ (e.g., equation (3.18)). For time-accurate 
calculations, < 7 y + i = a (a J + i) and can be expressed as 


°( z ) = -hAM-A* 2 ] > 0. 

For steady-state calculation and/or implicit methods (6 ^ 0 in (3.21)), 


(3.45e) 


< 7 ( 2 ) = -xj)(z) > 0. 


(3.45f) 


In other words, (3.45) is a first-order numerical flux with / replaced by / and the mean value charac- 
teristic speed CLj+k replaced by the modified characteristic speed a J+ i = a J + 1 + 7 J+ ±. Other more 
complicated forms of the g 3 function which include artificial compression can be found in Harten’s 
original paper. The current form (3.45) is quite diffusive, and a slight modification of this form without 
the use of additional artificial compression [ 6 ] by the author [43,44] will be discussed in section 3.5.2. 
To illustrate the difference in shock resolution between equations (3.45) and the form suggested by 
the author, numerical examples for one-dimensional shock-tube problems will be given in section IV. 

Roe-Sweby Second-order TVD Scheme : The scheme of Roe-Sweby starts with a first-order upwind 
scheme whose numerical flux is 


Vi = O I/ 3+1 + fj - s gnK + i)(/j+i - />■)]> 


(3.46) 


and adds a second-order correction term to hj+ 1 . The second-order numerical flux is of the form 

V i = h 3 + ± + -V[sgn(o i+ i) - Aa i+ i](/y+i - fj)}, (3.4 


where 


w j'-(-l+CT u j+c r 

T= A 3+i 


, o — sgn(a J+ i). 


Here 5(r) is a limiter and it can be 


S(r) — minmod(l,r), 

(3.47c) 


(3.47d) 

S(r) ~ max [ 0 , min( 2 r, l),min( 2 , r)] . 

(3.47e) 


The last limiter designed by Roe, nicknamed “superbee” [7], is the most compressive limiter among 
the three. It is especially designed for the computation of contact discontinuities. 

Osher and Chakravarthy TVD Scheme: Instead of using a MUSCL approach, Osher and Chakravarthy 
started with a one-parameter family of semi-discrete schemes with numerical flux 
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(3-48) 



~ h i+i - 






where = h(uy, uy+i) is some first-order upwind numerical flux. It can be the Engquist and Osher 

or Roe’s first-order upwind numerical flux. Here rj has the same meaning as before. 

The superscript 4- or - in ( ) denotes the flux difference across the wave with positive or negative 
wave speed. To obtain a higher-order TVD scheme, they modified the last four terms on the right-hand 
side by utilizing “flux” limiters as 


with 


Toe _ i 
3+4 _ h 3+\ 


(! *?) / a r-\ (* ^ ( \ 

—r~y a 3+v ) — T~y A 3+i f ) 


lL 7 J!i (A i+ */ + ) + iL r^(^-*/ + ), 


(3.49a) 


A J+ |/~ = minmod [A J+ |/ ,wA y+ i / ], (3.49b) 

A j+±f- = minmod[Ay + A/ _ ,u>Ay + |/~], (3.49c) 

A y+ i/+ = minmod[A y+ i/ + ,a;A J _i/ + ], (3.49d) 

Ay a /+ = minmod [Ay_ a / + , u> Ay + a / + ] , (3.49e) 


and 1 < u < I — 

— — 1-17 

One way to obtain a second-order time-discretization is to replace the forward Euler time- 
discretization by some linear multistep method or by the Runge-Kutta type of time-discretization. 
Note that the MUSCL way of obtaining higher-order in time is no longer valid in the Osher and 
Chakravarthy formulation. 


Due to the design principle of this scheme, for r) — — 1, the numerical flux h°2 L requires one more 
limiter than Harten’s method. For r\ ^ ~ 1, two more limiters than van Leer’s method and three more 
limiters than Harten’s method are required in addition to an extra step of getting higher-order time- 
discretization. The extra computations will become even more apparent as we extend these schemes 
to system cases. See reference [22] for details. 


3.5.2. Higher-Order Symmetric TVD Schemes 

The previous four methods are upwind methods. Next, the basic idea of second-order symmetric 
TVD schemes of Davis [8], Roe [9] and Yee [10] will be briefly described. Interested readers should 
consult the references cited for the actual construction of these methods. 

In 1984, Davis [8] expressed a particular form of Roe-Sweby’s second-order TVD scheme [7,41] as 
a sum of two terms. One term was the Lax-Wendroff method and the other term was an additional 
conservative dissipation term. He then simplified the scheme by eliminating the upwind weighting 
of the dissipation term and at the same time ensured that the simplified scheme still had the TVD 
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property. Shortly after that, Roe [9] reformulated Davis’s approach in a way that was easier to analyze 
and included a class of TVD schemes not observed by Davis. Subsequently, the author [45,46,10] 
generalized the Roe-Davis schemes to a one-parameter family of second-order explicit and implicit 
TVD schemes. The formulations of Roe-Davis can be considered as members of the explicit schemes. 
The main advantages of the author’s formulation are that stiff problems can be handled by using 
implicit methods and that steady-state solutions are independent of the time-step. 

A general discussion with extensive numerical examples can be found in a reference by Yee [10]. A 
careful examination of the modified-flux approach of Harten (later modified by Yee), and the symmetric 
TVD schemes of Davis, Roe, and Yee, reveal that these schemes have a very similar structure and 
can be expressed in the same general form. They are simpler to implement than the MUSCL or 
the Osher and Chakravarty schemes. Therefore, most of the numerical examples given later mainly 
employ these methods. Consider the general one-parameter family of explicit and implicit schemes of 
the form (3.21). The numerical flux for the second-order TVD schemes of Harten- Yee-Roe-Davis can 
be written as 


hj+k ~ 2 ^ + ( 3 - 50 ) 

The schemes only differ in the forms of the (f) function which are very similar to each other. 

Harten and Yee Upwind Scheme - Modified- Flux Approach : The <j> function of Harten ’s original 
modified-flux scheme was discussed earlier. That form of the numerical flux is quite diffusive. The 
author’s modification to equation (3.45a) is less diffusive and can be written as 


(3.51a) 

(3.51b) 
(3.51c) 

(3.51d) 

Here for an explicit method one sets ft — 1 , whereas for an implicit method or steady-state calculations, 
one sets ft — 0 and 9^0. This modified form is just a change in the definition of the original g 3 
function of Harten (equation (3.45b)) by removing the a J± i from (3.45b). In equation (3.51a), the 
is then incorporated as a factor of ( gj + <?y+i). 

One can generalize this method even further by including other limiters suggested by Roe and van 
Leer, such as 


‘t’j+i =^K + i)fc + 9 J+ i) - H a i+i +ij+j) A j+i. 


°( z ) = 

g, = minmod(A 3+ i,A y _i), 


'Ij+i =<7 K+i) 


gy+i-gy 




9i = S- max|o,min(2|A 3 . + i|,SA 3 ._i),min(|A 3+ i|,25A J _i)); 5 = sgn(A 3+ j ), (3.51e) 

(3.51f) 


9j = 




A,_ 1 + I A 


3 + 2 ^ 3- 2 


^ 3 + 5 ^ 3 - 

or even the simpler one, 


9 3 + 9 3 +i = minmod ( A i , A J+ 1 , A J + 1 ) , 


(3.51g) 
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with lj+i = 0. The minmod function of a list of arguments is equal to the smallest number in absolute 
value if the list of arguments is of the same sign, or is equal to zero if any arguments are of opposite 
sign. 


Yee-Roe-Davis Symmetric Scheme : For the Yee-Roe-Davis symmetric TVD schemes, the can be 
expressed in the form 


tj+i ~ [^P( a j+$) Qj+± + ^( a j+|)(^j+3 

with Qj + 1 chosen from 

= minmod(A J + i, Aj_±) + minmod(Ay + i, A J + |) - A j+ i, 
Qj+ 1 = minmod(Ay_i, Ay + i, Aj + |), 

Q j+ i =minmod[2A J _i,2A :)+ i,2A y+ |,-(A J _i + A i+§ )]. 


(3.52a) 


(3.52b) 

(3.52c) 

(3.52d) 


The coefficient of the first term on the right-hand side of <f>j+± is the second-order time-discretization, 
and the last term can be viewed as a numerical dissipation term. The parameter (3 has the same 
meaning as in (3.51). For the explicit method (0 = 0), if one sets (3 = 1 and Q to be the first limiter, 
it is the original explicit symmetric scheme of Davis. If one sets (3=1 and takes any of the three 
limiters, it is Roe’s TVD Lax-WendrofF scheme [9], Taking the implict scheme (0 ^ 0) with (3 = 0 and 
any of the limiters, it is the form that the author proposed. It is suitable for time-accurate as well as 
steady-state computations [46,10]. 

For analysis purposes it is sometimes convenient to let Qj+± 1 A i+ 1 and to express the 1 

function for the explicit second-order symmetric TVD schemes [9] as 


Vi = - [ A ( a ,-+±) S< ?j+i 

+ V>K+i)(i - <3,'+j)]Aj+i. 

(3.53a) 

- Vi 

r = . 

. .+ _ Vi 
Vj' 

(3.53b) 


The Q function can then be written as 


Q{r ,r + ) = minmod(l, r ) + minmod(l, r + ) - 1, (3.53c) 

<3(r~,r + ) = minmod(l, r“, r + ), (3.53d) 

<3(r - ,r + ) — minmod ^2, 2r~,2r + , ^(r“ + r+ )^ • (3.53e) 

The graphical representations of these three limiters for symmetric TVD schemes are shown in figure 
(3.3). Theoretically, one can design other limiters graphically by the aid of the sufficient condition. 


3.5.3. Global Order of Accuracy of a Second-Order TVD Scheme (Warming and Yee [47]) 

One of the drawbacks of higher-order TVD schemes is that they reduce to first-order at points 
of extrema. In the modified-flux approach, for example, the form of g 3 devised by Harten has the 
property of switching the second-order scheme into first-order at points of extrema (i.e., g 3 = 0 at 
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points of extrema). To see this, the behavior of the modified- flux approach around points of extrema. 
is examined by considering its application to data where 

Uj~i < Uj = Uy_|_i > • (3.54) 

In this case g 0 - — gj + 1 = 0 in (3.45c), and thus the numerical flux (3.45a) becomes identical to that 
of the original first-order-accurate scheme. Consequently, the truncation error of the second-order 
scheme (3.34) together with (3.45) deteriorates to 0((Ax) 2 ) at j and j + 1. This behavior is common 
to all TVD schemes, since this is one of the vehicles used to prevent spurious oscillations near a shock. 
Thus, a second-order TVD scheme must have a mechanism that switches itself into a first-order- 
accurate TVD scheme at points of extrema. Because of the above property, second-order-accurate 
TVD schemes are genuinely nonlinear; i.e., they are nonlinear even in the constant-coefficient case. 
Due to the uncertainty of the effect of the above property on the global order of accuracy, some 
numerical experiments were performed on the Harten scheme with an artifical compression [6] for 
Burgers’ equation 


du d . o . . 

Tt + Yx ( u/2 ) = °- ( 3 - 55a ) 

Here the flux function f(u) = u 2 / 2. Since the theory of TVD schemes is only developed for initial- 
value problems at this point, a periodic problem was considered to avoid extra complication. The 
initial condition is the same as shown in figures (3.1) and (3.2), namely 

u(x,0) = sin7ra;, 0 < x < 2. (3.55b) 

The local error of the computation at each grid point (jAz, nAt) is defined as 

ey = u” - u(jAx,nAt), (3.56) 

where u(jAx,nAt) is the exact solution of the differential equation (3.55). Here we assume that there 
is a fixed relation between Af and Ax. The global order of accuracy m is determined by 

||e|| = 0(Ax m ) (3.57) 


as the mesh is refined for some norm. 

To obtain the global order of accuracy numerically, the error at a fixed time was computed for 
a given mesh and repeated with increasingly finer meshes. Figure (3.4) shows the global order of 
accuracy of the second-order TVD scheme compared with the Lax-Wendroff method at time t — 0.2, 
when the solution is still smooth. The order of accuracy for the TVD scheme is 2 for the L x norm, 
around 3/2 for the L 2 norm, and 1 for the L norm. On the other hand, the order of accuracy for the 
Lax-Wendroff is 2 for all three norms. The main reason for the difference in the order of accuracy on 
the three norms for the TVD scheme is that the scheme automatically switches itself into first-order 
whenever extreme points are encountered. In this case there are two extreme points. 

Next, the global order of accuracy of the two methods was examined at time t = 1.0 when a shock 
has developed. Figure (3.5a) shows the order of accuracy of the TVD scheme at t = 1.0, which is 
identical to the one at time t = 0.2. But the order of accuracy for the Lax-Wendroff is drastically 
degraded. It is 1 for the L x norm, around J '2 for the L 2 norm, and 0 for the L norm. This is due to 
the inherent characteristic of the Lax-WenJ.roff method that causes this scheme to generate spurious 
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oscillations near the shock. Figures (3.1) and (3.2) show the numerical solution of the Lax-Wendroff 
method compared with the second-order explicit TVD method at t — 0.2 and t — 1.0 


3.6. Predictor-Corrector TVD Schemes to Include Source Terms 

All of the second-order explicit TVD schemes discussed so far are for homogeneous hyperbolic 
PDE’s. Consider a nonlinear nonhomogeneous hyperbolic conservation law 


du df(u) 
dt dx 


= s(u). 


(3.58) 


As noted at the end of section 3.3, the TVD property is only valid for the homogeneous part of 
equation (3.58). Certain types of source terms s(u) might preserve the original TVD property of 
the homogeneous part of (3.58), and others might not. However, disregarding the type of bounded 
source terms, one is not precluded from the use of TVD scheme when source terms are present, but 
precaution has to be taken in the procedure of including the source term. 

To include the source term efficiently, one can use (a) the method of lines with a linear multistep 
approach, (b) a two-step Lax-Wendroff type (e.g., explicit predictor-corrector MacCormack type), or 
(c) operator-splitting procedure (similar to the time-splitting procedure in multidimensional problems 
except the operator-splitting procedure is on the homogeneous part and the source terms). A discussion 
and derivation of the two-step Lax-Wendroff method can be found in reference [48]. In collaboration 
with Professor R. LeVeque (University of Washington, Seattle, Washington), research is underway 
to study the various ways of including stiff source terms in conjunction with the TVD property for 
time-accurate and steady-state hypersonic flows. 

For steady-state application, to avoid additional treatment of intermediate boundary condition and 
save storage, a straightforward way of extending the second-order explict TVD scheme to include source 
terms is to first rewrite the numerical flux without the source term in two parts: namely, a predictor- 
corrector Lax-Wendroff part and a conservative numerical dissipation part. One then includes the 
source term in the predictor-corrector step and considers the conservative numerical dissipation part 
as a second corrector step. Take for example the second-order explict symmetric TVD schemes ((3.50) 
together with (3.52)). The proposed predictor-corrector scheme can be written as 





+ Ate? 


(3.59a) 



+ U 1 ~ A 


/•( U 
U+l 



+ A( S « 1) } 


(3.59b) 


(3.59c) 

Here the superscripts “(1) and (2)” designate the values of the function evaluated at the intermediate 
solutions and u^ 2 K Also has a slightly different form than (3.52a), 


t /* +1 

3 


«?> + 


i#W-* 


■Ay+i = ~ "y+jK A y+§ - 

"y+i = Xa i+i- 
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(3.59d) 

(3.59e) 



where rp(z) is (3.18). The value + i can be any of the forms defined in equations (3.52). By defining 
a more complex , scheme (3.59) can be made upwind-weighted and would belong to the class of 
upwind schemes. The derivation is straightforward and will not be given here. 

One can see that the formulation of this scheme is broken into two parts, namely, the predictor- 
corrector step of the MacCormack explicit scheme, and an appropriate conservative dissipation term. 
Here the predictor-corrector scheme is TVD in the sense of a constant-coefficient homogeneous case 
(s = 0) and with 4>j±± evaluated at u n instead of t/ 1 ) . For the general nonlinear case, it appears 

to be difficult to prove that this predictor-corrector scheme is TVD, but numerical experiments for 
one and higher-dimensional nonlinear homogeneous hyperbolic conservation laws show that (3.59) 
has the TVD-type properties. Other equivalent predictor-corrector forms can also be used. This 
predictor-corrector TVD method is sometimes referred as the “TVD MacCormack” scheme. It is a 
slight modification of Roe’s one-step TVD Lax-Wendroff scheme. If one sets Q to be equation (3.52b) 
and tp(z) = \z\, the scheme is the same as described in Davis [8] and Kwong [49]. The reason for 
choosing the predictor-corrector step instead of the one-step Lax-Wendroff formulation is that the 
predictor-corrector form provides a natural and efficient inclusion of the source terms especially for 
multidimensional problems [48]. 


3.7. Semi-Implicit TVD Schemes for Problems Containing Stiff Source Terms 

The explicit TVD scheme (3.59) can be used for either time-accurate or steady-state calculations. 
It is second-order accurate in time and space. However, for time-accurate calculations, (3.59) is TVD 
in the sense of the constant-coefficient homogeneous case and with <j> ■ . i evaluated at tx n instead of 

u^ 1 '. Moreover, if the source term is stiff, the restriction in the time-step due to stability requirements 
is prohibitively small, and (3.59) is not practical, especially for steady-state applications. In this 
section, a semi-implicit method is proposed for steady-state computations. Another alternative is a 
fully implicit method. The basic implicit scheme and the related difficulty in extending the implicit 
method to higher dimensions with stiff source terms will be discussed in later sections. 

The idea of treating the stiff term implicitly and the non-stiff term explicitly is a common procedure 
in numerical methods for stiff ordinary differential equations. The semi-implicit treatment for PDE’s 
with stiff source terms in conjunction with classical shock-capturing method is also a common proce- 
dure; see for example, reference [50]. What is proposed here is to replace the classical shock-capturing 
methods with a modern shock-capturing method. If one follows the idea of Bussing and Murman [50] 
in treating the source term implicitly, a semi-implicit predictor-corrector TVD scheme can easily be 
obtained. The basic idea is to treat the source term implicitly and the homogeneous part of the PDE 
with a predictor-corrector TVD scheme. For “extremely” stiff source terms, it is advisable to solve 
the resulting nonlinear system iteratively. However, for a “moderatelly” stiff source term, in order 
to avoid solving nonlinear equations iteratively, the Taylor expansion of the source term at time-level 
n + 1 is truncated to first-order as 


^ +1sSS " + (£)"K +I - u ?)- (3-60) 

The scheme can be written as a one-parameter family of time-differencing schemes for the source term; 
i.e., the following formulation includes scheme (3.59). The proposed scheme is 

d l Au l' ] = (/" -/"-■) + Ats 7> (3.61a) 
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(3.61b) 


'= ( x -*"!;)• 

Au< 2) = - A - j f> + A«4 IJ ), (3.61 c) 

U "+ 1 = U < 2 >+ , (3.61d) 

with + Uy and u,p = + Uy 1 ^. Here, d is assumed to be nonzero; i.e., only the 

type of source terms such that d is invertible at each grid point are permissible. The parameter 9 is 
in the range 0 < 9 < 1. For 0/0, the source term is treated implicitly. If 9 = 1, the time-differencing 
for the source term is first-order, and (3.61) is best suited for steady-state calculations. Note that the 
order of time-accuracy, which is determined by the parameter 9, has a different meaning than for the 
9 appearing in the implicit method (3.21). To obtain a second-order time-discretization, one can set 
9 — 1/2, and (3.61c) is replaced by 

■ 

d? Auf 1 = -A /]*>, - fj 1] + At*?. (3.61e) 

and 

„(. 2 > = U? + ^(Au* 11 + A«' J) ) (3.611) 

Equation (3.61e) is very similar to (3.61c) except d and s are evaluated at u n instead of t/ 1 ), and 
uP is (3.61f). By doing this, scheme (3.61,a,b,e,f) is second-order in time and space (R. LeVeque, 
University of Washington, Seattle, Washington, private communication). 

3.8. Linearized Form of the Implicit TVD Schemes 

All of the TVD schemes discussed above are nonlinear schemes in the sense that the final algorithm 
is nonlinear even for the constant-coefficient case. For implicit TVD schemes (9 / 0 in equation 
(3.21)), the value of u n+1 is obtained as the solution of a system of nonlinear algebraic equations. To 
solve this set of nonlinear equations noniteratively, a linearized version of these nonlinear equations is 
considered. For the non-MUSCL formulation, linearized forms can easily be obtained. For illustration 
purposes, only the linearized form of implicit symmetric TVD schemes will be discussed. The same 
idea can be used for the implicit upwind TVD scheme (3.51). A detailed derivation can be found in 
Yee [43]. Also, unlike the Lax- Wendroff- type scheme, it is more straightforward to include the source 
terms for the implicit scheme (3.21). See section 6.4 for a discussion. 

3.8.1. Linearized Version for Const ant- Co efficient Equations 

For the linear scalar hyperbolic PDE (2.1), the numerical flux for the symmetric TVD scheme 
together with (3.53) can be written as 

| [ a («j+i + u j) - l°l( 1_< 2,+i)A 3+ i]. (3.62) 

Substituting (3.62) in (3.21), one obtains 
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u n+1 + — 

J 2 

_ A0 
“ 2 


au 


3 + 1 ~ W(1 - Q j + i)A j + lU 


n + 1 


1 n+1 


au ,- 1 ~ \a\{l-Q j _i)A ;j _i 


RHS of (3.21). 


(3.63) 


Here “RHS of (3.21)” means the right-hand side of equation (3.21) with h n . defined in (3.62). Locally 

3~T 2 

linearizing the coefficients of (A J± i«) n+1 in (3.63) by dropping the time-index from (n-f 1) to n , one 
obtains 


u" +1 + — 

5 2 


au 


n+1 


M«?-i -|a|(l-Q? +i )A 3 . + i u 


n+l 


J + • 


+ |a|(l - i )A •_ i u n+1 l = RHS of (3.21). 


Letting d 3 — u^ +1 - u” (the “delta” notation), equation (3.64) can be written as 


eidj-i + e 2 dj + e 3 d y+1 = -A (^” + i - 


where 


(3.64) 


(3.65a) 


ei 




e 3 


A0 

2 


X9 

1 + y 


A/9 

2 


a- |a|(l-Q y _i) 

fl|(l-Q J -_i) + |«|(l-Q i+ i) 

r 

a - |a|(l ~Q 3 +l) 


(3.65b) 

(3.65c) 

(3.65d) 


The linearized form (3.65) is a spatially five-point scheme and yet it is a tridiagonal system of linear 
equations. This is because at the (n+ l)th time-level only three points are involved; i.e., u n j' 1 1 ,u n+1 , 

l < 3 3 

and Although the coefficients involve five points, they are at the nth time-level. 

The form of (3.65) is the same as (3.64) except the time-index for the Q ± i and ^ , is dropped from 

(n + 1) to n for the implicit operator. One would expect that the linearized form (3.65) is still TVD. 
Numerical studies on one- and two-dimensional gas-dynamics problems supported this hypothesis. It 
was found in reference [51] that when time-accurate TVD schemes are used as a relaxation method 
for steady-state calculations, the convergence rate is degraded if limiters are present in the implicit 
operator. Therefore, for steady-state applications, one might want to use the linearized form obtained 
by setting Qj±i — 0 in (3.65); i.e., by redefining the coefficients in (3.65) as 


^2 = 1 + A#(|a|) , 

Xe ( Ml 

e s = y(«“ M)- 
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(3.66a) 

(3.66b) 

(3.66c) 



Scheme (3.65a) together with (3.66) is spatially first-order accurate for the implicit operator and 
spatially second-order accurate for the explicit operator. 

3.8.2. Linearized Version for Nonlinear Equations 

For the nonlinear case, the situation is slightly more complicated since the characteristic speed 
df fdu is no longer a constant. For the symmetric TVD scheme, after substituting (3.53) in (3.21), 
one obtains 


n+l 



fj + 1 

Si - 1 




n+l 


n+l 


= RHS of (3.21). 


(3.67) 


Unlike the constant-coefficient case, one also has to linearize ^(°"±1)> an ^ Following 

the same procedure as in [10,43], two linearized versions of (3.67) are considered. 

Linearized Nonconservative Implicit Form : Adding and substracting /™ +1 on the left-hand-side of 
(3.67) and using the relation (3.14), one can express (3.67) as 


xe 
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n+l 

l i+i 






‘J + i 


l u 


n + l 


J 2 •* 2 J 2 


A, iu 

J 2 


n+l 


= RHS of (3.21). 


(3.68) 


Rewriting (3.68) in the same form as (3.31) and dropping the time-index of the coefficients of 
u ” - *" 1 from in + 1 ) to n , one obtains 

e\dj-\ + 62 dj + e% dy+i = — h™_ (3.69a) 


where 


and 


e x = A OB", 

(3.69b) 

ej = l-A#(s- + B + ), 

(3.69c) 

e 3 = \6B + , 

(3.69d) 

±a j±i - V>(a 3 ±i)(l-Q 3 ±i) 

n 

(3.69e) 


Equation (3.69) is again a five-point scheme, and yet the coefficient matrix associated with the dj’s is 
tridiagonal. With this linearization, the method is no longer conservative. Therefore, (3.69) is more 
applicable to steady-state calculations. A spatially first-order- accurate implicit operator similar to 



(3.69e) can be obtained for (3.69) by setting B ± = |[±o i=t i - rp(±a j± i)] n . Since the limiter does 
not appear on the left-hand side, improvement in efficiency over (3.66) might be possible [43,24]. This 
reduced form is especially useful for multidimensional, nonlinear, hyperbolic conservation laws. 


Linearized Conservative Implicit Form : One can obtain a linearized conservative implicit form by 
using a local Taylor expansion about u n and expressing f n+1 - f n in the form 

fr - fj = <*?K +1 - w ?) + <>(At 2 ), (3.70) 

where a" = (df /du)”. Applying the first-order approximation of (3.70) and locally linearizing the 
coefficients of (A Jrt iti) n+1 in (3.67) by dropping the time-index from (n + 1) to n, one obtains 




n + 1 


+ 


XO r 
2 L 


+V>(a"_.)(l - A y _iu n+l ] = RHS of (3.21). 


n+l 


Letting dj = u^ +1 — u”, equation (3.71) can be written as 


e \ dj-i + e^dj + e s dj+i — -X (h“ + i - h ”_ A ), 


(3.71) 


(3.72a) 


where 


(3.72b) 
(3.72c) 
(3.72d) 

The linearized form (3.72) is conservative and is a spatially five-point scheme with a tridiagonal system 
of linear equations. Scheme (3.72) is applicable to transient as well as steady-state calculations. As 
of this writing, the conservative linearized form (3.72) has not been proven to be TVD. Yet numerical 
study shows that for moderate CFL number, (3.72) produces high-resolution shocks and nonoscillatory 
solutions. 


\e 


e x = 


a 3-i - V>K_i)(i -Qi-k) 


e 2 = 1 + 


xe 


^( fl j-i)(l - Qj-h) + - Q,-+i) 


: xe 


es 


a j+ 1 ^( a j+f)(^ Qj+0 


For steady-state applications, one can use a spatially first-order implicit operator for (3.72) by 
simply setting all the Q J± i = 0; i.e., redefining (3.72b)-(3.72d) as 


ei = 


XO 


e 2 = 1 + 


xe 




V > (a i _i) + 0(a y+ i) 


e 3 = 


XO 


a 3 + 1 ~ 


(3.73a) 

(3.73b) 

(3.73c) 


Numerical experiments with two-dimensional steady-state airfoil calculations show that this form 
(alternating direction implicit (ADI) version) is the most efficient (in terms of CPU time) among 
the various proposed linearized methods for the case of 6 — 1. No comparison has been made for 
time-accurate calculations or for any othe values of 0 < e < 1. 
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IV. Extension of Nonlinear Scalar TVD Schemes to 1-D Nonlinear Systems 

Before getting into a detailed discussion of this section, it is important to emphasis that first-order 
upwind schemes for one-dimensional nonlinear systems encountered in the literature differ mainly 
by their so called “Riemann solver”. For constant-coefficient systems, they all reduce to the CIR 
(Courant- Isaacson- Rees) method [36] . There exists three popular ways of extending scalar schemes to 
nonlinear systems (hereafter, referred to as Riemann solvers): the exact Riemann solvers [28,52], the 
approximate Riemann solvers [53,54], and the flux-vector splitting techniques [55,56]. In this paper, an 
approximate Riemann solver of Roe and the flux-vector splittings of Steger and Warming and of van 
Leer for perfect gases are considered. The generalized Roe’s approximate Riemann solver of Vinokur 
[57], and the generalized flux-vector splittings of Vinokur and Montagne [58] for real gases are also 
considered. 

Recent developments of higher-order modern shock-capturing methods stress on the nonlinear scalar 
case. The extension of higher-order modern shock-capturing methods relies heavily on these Riemann 
solvers. Furthermore, since the extension is not unique, it depends on the form of the nonlinear schemes 
that one started with [22]. One form might be simpler than the other for its system counterpart. The 
situation arises even when one starts with a scheme with two different representations for the scalar case 
(see reference [22] for details). For example, the Osher and Chakravarthy scheme is more complicated 
and more expensive in the system cases than other schemes under discussion. Moreover, a comparison 
among the numerical results for schemes (3.39), (3.47), (3.49), (3.51), and (3.52) does not indicate 
any advantage of the more complicated schemes over the simpler ones. Based on this fact, numerical 
results presented here reflect the author’s personal experiences and preference for certain schemes. 
No attempt has been made to present a unified comparison. Also, no effort has been made to collect 
numerical results from investigators in related fields to illustrate the performance of similar schemes. 
Readers are encouraged to study the related theory and numerical results of references [6,59-68]. 

Since the current discussion is on conservative shock-capturing finite-difference methods, noncon- 
servative schemes such as the A-scheme of Moretti will not be discussed, although recent results of 
Moretti et al. show that the A-scheme together with a shock-fitting procedure suggests an efficient 
alternative to shock-capturing methods. More detailed study and numerical tests are needed in this 
direction in the near future. Also, the FCT scheme of Boris and Book [38] for nonlinear systems does 
not make use of any type of Riemann solvers and thus will not fall into the category of the current 
discussion. Historically, Boris and Book were two of the pioneers in introducing the concept of flux 
limiters for nonlinear scalar hyperbolic conservation laws. However, for nonlinear systems they applied 
flux limiters to the individual flux functions. In the author’s opinion this is less effective than the use 
of Riemann solvers. A discussion is included below. Also the popular and efficient shock-capturing 
method of Jameson et al. will not be discussed here. The scheme of Jameson et al. [18] seems to fall in 
between the classical and modern shock-capturing schemes. Their scheme, originally designed for the 
Euler equations, consists of two adjustable parameters to control the amount of numerical dissipation, 
but does not use any Riemann solver. The relative advantages and disadvantages of the schemes 
of Moretti, Boris and Book, Jameson et al., and TVD schemes for weak and moderate shock-wave 
calculations are not clear at this point. However, for strong shock waves, especially in the hypersonic 
regime, TVD schemes in conjunction with the appropriate Riemann solvers are conjectured to perform 
better than the other aforementioned approaches. 

4.1. Methods of Extension (Riemann Solvers) 

The following is a brief review of the methods of extending nonlinear scalar difference schemes to 
nonlinear systems of hyperbolic conservation laws. The objective is to give a flavor of the developments 
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and therefore many of the details are left out. First the various methods of extension to systems and 
the original use of these methods in conjunction with the first-order upwind finite-difference methods 
will be discussed. All of the original methods of extension to systems were developed for perfect or 
thermally perfect gases and the resulting algorithms for the gas-dynamics equations were first-order 
accurate (except in flux- vector splitting approaches). In the subsequent section, generalization of these 
methods to real gases will be described in conjunction with numerical schemes that are higher than 
first-order. 


The conservation laws for the one-dimensional Euler 


equations can be written in the form 


dU dF(U) 
dt dx 


= 0 , 


where the column vectors U and F(U) take the form 


(4.1a) 



p 


pu 

u = 

rrt 

, F = 

mu + p 


e 


eu + pu 


(4.1b) 


Here p is the density, rrt — pu is the momentum per unit volume, p is the pressure, e = p(c + |u 2 ) is 
the total internal energy per unit volume, and e is the specific internal energy. 


Many approximate Riemann solvers make use of the eigenvalues and eigenvectors of the Jacobian 
matrix A — dF/dU. For a general gas, one therefore requires the thermodynamics derivatives of p. 
In terms of the internal energy per unit volume 7 — pe , the thermodynamic derivatives can be defined 
as 


Here the subscript 7 
( p ) constant. If h 
relation 


dp 

dp 


dp 

dc 


(4.2) 


( p ) means the partial derivatives of p with respect to p (?) by holding 7 
e + p/p is the specific enthalpy, one can obtain for the speed of sound c the 


c 2 = x + K h- 


For a perfect gas, x — 0, and k, = (7 - 1). 

The Jacobian matrix A takes the form 

0 1 0 
X — (2 — k)u 2 / 2 (2 - /c)u k 

(x + /cu 2 /2 — H)u H — ku 2 (1 + k)u 

where H = h + u 2 /2 is the total enthalpy. The three eigenvalues of A are 

a 1 = u — c, a 2 — u, and a 3 — u + c. 
The corresponding right-eigenvector matrix is 




1 

u — c 
H-uc 


1 

u T c 
K H + uc 


(4.3) 


(4.4) 


(4.5) 


(4.6) 
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while its inverse can be written as 


where b i 
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-b 2 
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- i(«*« - i) 

Li 
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(4.7) 


In order to relate the variables p and c to the independent variables p and e, it will also be convenient 
to introduce the nondimensional thermodynamic variables 


1 = 




(4.8) 


For a perfect gas, these two parameters are constant and equal to each other; for an equilibrium real 
gas, they are both arbitrary functions of p and e. 


Again, many existing conservative schemes for the system (4.1) use forward Euler time-discretization 
and the scheme can be written as 


[/;+' = u* _ a [?; + i -F;_ 4 ]. (4.9) 

The vectors Fj ± i are numerical flux vectors corresponding to h JZ t i for the scalar case. Note that the 
numerical flux function F J+ 1 should not be confused with the flux function Fj. 

CIR Method : The earliest method for gas-dynamic equations in characteristic form was proposed 
by Courant et al. [36]. Their procedure, sometimes called the CIR method, is to trace back from 
(j Ax, (n -f 1)A t) all three characteristic paths. Since the problem is nonlinear, the directions of these 
path are not known exactly, but to a first approximation they can be taken equal to their known 
directions at ( jAx,nAt ). Then each characteristic equation is solved using interpolated data at time 
rcAt, in the interval to the left of j for characteristics with positive speed, and in the interval to right 
of j for characteristics with negative speed. The resulting method is only first-order accurate. It has 
a principal drawback in that the scheme cannot convey a shock wave with the proper speed because 
it is not a conservative scheme. This method was later rediscovered by Chakravarthy et al. [69] and 
was renamed the split-coefficient matrix method. 

Exact Riemann Solver : Godunov [28] was the first to develop the idea of advancing the solution to 
the next time-level by solving a set of Riemann problems. Recall that the Riemann problem for any 
system of conservation laws arises if initial data are prescribed as two semi-infinite states ( U = U L 
for x < 0 , U = U R for x > 0). The solution then consists of centered waves. For the one-dimensional 
Euler equations, they consist of three waves (u,u ± c); the inner one is a contact discontinuity, and 
the outer ones may be shock waves or rarefaction fans. The exact solution of this problem involves 
only algebraic equations. See reference [19,28] for details. 

Let U n j be the average state over ((y ± |)Az) at (nAt). The way Godunov used the Riemann 
solvers was to replace the data by an approximate distribution in which the state inside each interval 
is uniform and equal to U™ . For each interface (j + |)Aa;, one can so ^ ve the Riemann problem 

with U L — Uj and U R = Uj + This gives an exact solution u U n +±” to the approximate problem, 

3 ' 2 

assuming At is small enough that the waves from neighboring interfaces do not intersect. The solution 
at (n+ l)At can again be approximated by a piecewise uniform distribution, and then the process can 
be repeated. For the Godunov method, the numerical flux is 
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(4.10) 



= *Kt 


)• 


More recently Ben-Artzi and Falcovitz [70] generalized the exact Riemann solver of Godunov to be 
second-order accurate. Their theory is too complicated to be summarized here. The versatility of 
their method remains to be shown. 

Another closely related method, devised by Glimm [71], modified by Chorin [72], and further im- 
proved by Colella [52], is the random-choice method. It represents the staggered grid solution by 
randomly sampling the Riemann solutions generated at the previous time-level. Recently, Toro [73], 
Roe and Toro [74], and Toro and Clarke [75] have extended the random-choice method to be second- 
order accurate. Their preliminary result is very encouraging for utilizing this type of scheme in 
combustion-related flows. 

Approximate Riemann Solvers : Since the Riemann problems arising in Godunov’s method relate only 
to an approximation of the data, one might reasonably be satisfied with approximate solutions of the 
Riemann problem if these solutions still describe the important nonlinear behavior. Roe [53] , Harten 
et al. [54], and Osher and Solomon [76] proposed methods for finding such solutions. The method of 
Osher and Solomon is quite complicated to explain, and uses similar arithmetic operation (and is not 
as exact) as the Godunov method. The method of Harten et al. does not retain all the information 
about every wave. Therefore, only Roe’s approximate Riemann solver is described here. Interested 
readers should refer to the original references for details. Roe’s approximate Riemann solver is a linear 
wave decomposition in which he required that there exists an average state U which is a nonlinear 
function of the left and right states U L,R satisfying 

(1) F(U r ) - F(U l ) = A(U r ,U l )(U r - U L ) = A{U)(U r - U L ) 

(2) A(U R , U L ) has real eigenvalues and a complete set of eigenvectors 

(3) A(U,U) = A(U) 

(4) U R - U L = R{U r , U L )a = R{U)a 

where the fcth element of a is the strength of the kth characteristic wave (or the jump in the char- 
acteristic variables). The main feature of the method that makes it valuable for nonlinear systems 
is that it returns the exact solution whenever U L and U R lie on opposite sides of a shock wave or 
a contact discontinuity. For the one-dimensional Euler equations for a perfect gas, the Roe averaged 
state can readily be obtained as 


with 


u L + Du R 

u = , 

1 + D ’ 

(4.11a) 

— h l + dh r 
H - 1 + D ’ 

(4.11b) 

= ( 7 - 1) [ff - iu 2 | , 

(4.11c) 





(4. lid) 
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To make use of his approximate Riemann solver, Roe applied his average state (4.11) in conjunction 
with his first-order upwind scheme for the nonlinear scalar case to arrive at a first-order scheme for 
the nonlinear system case. The numerical flux is of the form 

where Rj+± = R{U J + ±) and 


(4.12b) 


«,-+i = r ;+a u ^ - Uj). 

The subscript “j + denotes the average state between j and j + 1 using relation (4.11). 

Flux- Vector Splittings : The simplest way of introducing upwinding into systems of hyperbolic conserva- 
tion laws seems to be based on the representation of the flux vector F(U) as the sum of F~(U) + F + (U) 
such that one can apply forward- and backward-differencing on the Jacobian matrices and 
respectively. This would be equivalent to using an approximate Riemann solution in which the nu- 
merical flux F J+ a is of the form 

F, + , = F+^ +i ) + F -{Uf + ,). (4.13) 

This numerical flux amounts to requiring that the Jacobian matrices ^ — and — ^ have no positive and 
no negative eigenvalues, respectively. A popular way is to split the flux according to the characteristic 
speeds (u, u ± c). 

This idea, assumming that the flux has the homogeneous property of degree one (for thermally 
perfect gases), seems to have been first used in the context of astrophysical gas dynamics [77], and 
to have been rediscovered with a fuller mathematical development by Steger and Warming [55]. Note 
that Steger and Warming make use of the flux-vector splitting for a second-order upwind scheme in 
a non-MUSCL way and without the use of limiters (non-TVD method). Their final second-order 
upwind scheme cannot be represented in the form (4.13). Mulder and van Leer [78] and Anderson 
et al. [61] introduced the MUSCL approach with limiters into the Steger and Warming flux- vector 
splitting (TVD method) , and its numerical fluxes have the form (4.13). The MUSCL approach with 
limiters (TVD) formulation can produce a better shock resolution than the non-MUSCL and non-TVD 
approaches. 

In 1982, van Leer [56] devised an alternative splitting for a perfect gas such that there are noticeably 
better results around sonic points and sharper shock transitions than can be obtained with the Steger 
and Warming splitting. In the next few sections, two generalizations of these flux- vector splittings and 
a generalization of Roe’s average to real gases will be described. In these formulations, the perfect-gas 
version is included as a particular case. 


(4.12c) 


A m = 


o 

0 


a 2 : 
3 + h 




(4.12a) 


4.2. Description of the Riemann Solvers for Real Gases 

Extensions of the exact Riemann solver of Godunov to real gases have been obtained by Colella and 
Glaz [79], Dukowicz [80], and Ben-Artzi and Falcovitz [70]. The derivations are quite involved and 
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interested readers should refer to their papers for details. The recent generalization of Roe’s average 
by Vinokur [57], and generalization of Steger and Warming and van Leer flux- vector splittings by 
Vinokur and Montagne [58] will be discussed in this section. These generalizations are also simpler to 
use than the exact Riemann solvers of references [70,79,80]. 

The following three subsections (4.2.1 - 4.2.3) were written by M. Vinokur of Sterling Software, Palo 
Alto, California. They are summaries of his proposed manuscripts that are in preparation. 


4.2.1. An Approximate Riemann Solver (Generalized Roe Average [57]) 

Among the various approximate Riemann solvers [53,54] for a perfect gas, the most common one 
uses the Roe average [53] because of its simplicity and its ability to satisfy the jump conditions. The 
only drawback is that Roe’s “property U conditions” cannot be uniquely satisfied for nonperfect gases. 
Various alternatives were proposed for real gases [30,81,54], All the numerical examples illustrated 
later employ the arithmetic mean average similar to the ones suggested by Huang [30] and Carofano 
[81], and a form derived by Vinokur [57], The results of the numerical examples indicate that there is no 
drastic difference between the three averages for the one-dimensional shock-tube problem with various 
ranges of density, pressure and Mach number. The exact formulae of Vinokur will be given below, 
whereas the others can be found in the appropriate references. The reason for choosing Vinokur’s 
formulation is that his derivation is a more systematic approach than those of Huang and Carofano. 
The use of the approximate Riemann solver in conjunction with the numerical schemes will be discussed 
in sections 4.4 and 4.5. 

The flux at a point separating two states U 1 ' and U R is based on the eigenvalues and eigenvectors 
of some average A. The optimum choice for A is one satisfying 

A F = AAU , (4.14) 

where A(-) — (-) R — (-) L . This^ choice of A will capture discontinuities exactly. One way of obtaining 
A is to seek an average state U , such that 

A = A{U). (4.15) 

Such a state is known as a Roe-averaged state. Expressions for a perfect gas were first devised by Roe 
[53] and are given by equation (4.11). 

The entries in A depend explicitly on the thermodynamic variables h, x, and k. Since the density 
is not explicitly required, one would expect the Roe-averaged state to depend on p L and p R through 
their ratio only. It is therefore convenient to define the parameter D = \fp R / p L . We first examine 
the second component of equation (4.14). The average velocity u must be a linear combination of 
u L and u R . Recalling that u L and u R can be independently prescribed, we can readily establish the 
same u as in equation (4.11a) for a perfect gas. This definition will satisfy all the terms involving the 
velocity. Note that u always lies between and u R . The remaining terms in the equation result in 
the condition 


X A p + kAc = A p. (4.16) 

This last condition is automatically satisfied for a perfect gas. 

In order to satisfy the third component of equation (4.14) we also require H to have the same form 
as the perfect-gas version (equation (4.11b)). Using the definition of H, equations (4.11a) and (4.11b) 
can be combined to define the Roe-averaged specific enthalpy as 
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D 


(4.17) 


h = 


h L + Dh R 
1 + D + 


2(1 + D) 2 


(Au) 2 . 


Note that h could lie outside the range given by h L and h R if A u is sufficiently large. The Roe-averaged 
sound speed is given by (4.3) as 

c 2 = x + (4.18) 

For a perfect gas, equations (4.11a), (4.11b), (4.17) and (4.18) are sufficient to define uniquely a\ R , 
and R 1 , since x — 0 and « are given constants. 

For a real gas, equation (4.16) provides only one relation for the variables x and k. Since the 
equation of state can be non-con vex, a universally valid unique solution may not be possible. In order 
to gain some insight, we consider the special case in which the states R and L are precisely those that 
satisfy the jump conditions across a discontinuity (either a contact discontinuity or a shock wave). 
Then equations (4.11a) through (4.18) are consistent with the exact Riemann solver, even though x 
and k are not uniquely defined. For a shock wave one obtains 


and 


- h L + D 2 h R 

1 + D 2 


(4.19) 


-2 _ A P 
A p' 

One can show that no further relations are required in Roe’s approximate Riemann solver discussed 
in section 4.3. For the special case of a thermally perfect gas, c 2 is a function of h only, and one can 
readily show that the values of h and c 2 given by equations (4.19) and (4.20) can only satisfy this gas 
law if the function is linear. But this is precisely the definition of a perfect gas. 

The above analysis makes it clear that for a real gas the values of xjand #c must be defined in terms 
of the thermodynamic states R and L, and not in terms of the state h. Due to the nonconvex nature 
of a real-gas equation of state, the values of x and k at states R and L, or some average of the two, 
will not satisfy equation (4.16) in general. One way to obtain unique values of x and /c is to project 
the average state given by R and L onto the straight line defined by equation (4.16). Since the value of 
Ae depends on the arbitrary constant in the definition of e, the resulting value of k will depend on the 
choice of this constant. To overcome this arbitrariness, one divides equation (4.16) by k. The straight 
line for the variables 1 /k and x/« has a slope given by A p and A p, both of which are uniquely defined 
by states R and L. Since A p and A p are not dimensionally consistent, one must further introduce a 
scale factor proportional to their ratio in order to have the final answer independent of the choice of 
dimensional units. A natural candidate is the square of an average sound speed. In terms of arithmetic 
averages, and the scale factor c 2 = [(c L ) 2 + (c jR ) 2 ]/2, one obtains the expressions 

= = {£[(£ + S> Ap+ A7A *} / (V + 2W) (4.21) 

and 

I = {|[ ( S + $ )Ap+ + ^ A?Ap } / < Ap2 +2<Ap2 )- ( 4 - 22 > 
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Other expressions can be obtained by using different scalings and taking different types of averages. 
The optimum choice which would be valid for a wide range of conditions will require further research 
and numerical experiments. 


4.2.2 Generalized Steger and Warming Flux- Vector Splitting [58] 

For subsonic flow, the eigenvalues of A are of mixed sign. In flux-splitting methods, the flux F is 
divided into several parts, each of which has a Jacobian matrix whose eigenvalues are all of the same 
sign. The present class of flux splitting makes direct use of the eigenvalues of A , and is an extension 
of the work of Sanders and Prendergast [77] and Montagne [82]. The basic idea is the observation 
that the eigenvalues are actually three velocities. We can associate with each eigenvalue a l , l = 1,2,3, 
a stream with velocity a z , and some unknown density p l and specific internal energy € l . Each stream 
can then be characterized by the column vector 


(4.23) 


U l = 


m 


where m l — p 1 a 1 and e l = p l [t l + |(a z ) 2 ]. The flux due to each stream is assumed to be convective 
only, namely, 

F l = a l U l . 


The six unknowns are determined from the conditions 

3 


u = '^2 u ‘ an<J f =^ 2 f ‘. 


(4.24) 


(4.25) 


1 = 1 


1 = 1 


Since the second component of U and the first component of F give the same equation, we are left 
with one degree of freedom. From the first two components of U and F one readily obtains 


P 1 = P 3 = -f- and p 2 - p( 1 - -). 
27 7 


(4.26a) 


The third components of U and F result in the relations 


e 1 ^ 3 


and = 


(4.26b) 


1=1 


It is convenient to define the parameter a as the fraction of the internal energy per unit volume 
contained in stream 2. Therefore 


22 ( 3 “ 7) 

H = ape ~ 


and p 1 e 1 = p 3 e 3 = (1 - a)pe 


( 3 - 7 ) 


(4.26c) 


The final expressions for the F l can be written in the form 


F 2 = 


27 


2(7 - l)a 2 

2 ( 7 -l)(a 2 ) 2 

( 7 -l)(a 2 ) 3 + a^la 2 c 2 


(4.27a) 
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while 


^ + d-«) 


2 (^- 1 ) 


for / = 1 or 3. Note that a 2 here means the second component of the eigenvalue of A whereas (a*) 2 
means (a z )(a z ). We thus have a one-parameter family of flux splittings, where a can be an arbitrary 
parameter. The total flux can be split according to the sign of the eigenvalues. For -c < u < 0 we 
therefore have 


F+ = 


while for 0 < u < c, 


F 2 + F e 


= F 1 + F 2 , 


F~ = F 1 . 


(4.28a) 


(4.28b) 


For a perfect gas one can show that when -c < u < 0 the determinant of the split-flux Jacobian 
A+ = dF+/dU is 


det(A + ) — 


(3-7 )(c + u) 3 (l-a) 


(4.29a) 


The determinant is the product of the three eigenvalues, and a necessary condition for A + to have 
eigenvalues that are all positive is that det(A + ) > 0. It follows from equation (4.29a) that we must 
take 7 < 3. For the region 0 < u < c, one can show that the minimum of det(A + ) occurs when u 
approaches zero. This minimum value is 


3 

det(A + ) = ^3 5 — 37 — (3 — 7) 2 a . 


(4.29b) 


The largest value of this determinant is given by a — 0. For this case the determinant will be positive, 
provided that 7 < 5/3. By examining the other coefficients in the characteristic polynomial, one 
can show that all the eigenvalues are positive under these conditions. One can further show that in 
general the three eigenvalues of A + are discontinuous at u — — c, 0, and c. For a = 0, the two largest 
eigenvalues are continuous at u = c. Based on these considerations for a perfect gas, the value a — 0 
is the recommended value for a real gas. For a perfect gas it reduces precisely to the Steger- Warming 
flux splitting [55]. 


4.2.3 Generalized van Leer Flux- Vector Splitting [58] 

In a different approach, van Leer [56] constructed a flux splitting for a perfect gas in terms of low- 
order polynomials of u which gives continuous eigenvalues at u = 0 and u — ±c. The splitting also has 
the desirable property that one of the eigenvalues of the split-flux Jacobians is identically zero. This 
results in a sharper capture of transonic shocks. An extension of van Leer’s splitting for a real gas 
was derived by Montagne [83], but it is not internally consistent. The present formulas represent the 
most general, consistent extension of van Leer’s splitting for a real gas. Due to the arbitrary nature 
of a real gas law, the condition of one eigenvalue being identically zero cannot be satisfied exactly. 

For |u| < c, the continuity requirements necessitate a factor (u ± c) 2 in the formulas for F±. The 
expressions for the first two components of F ± that are given by the lowest order polynomials in u 
are readily found to be 


F ± = ±- 

x mass 


(u ± c)‘ 


(4.30a) 
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and 


74 c 


(it ± c) 2 [2c ± (7 — l)u . 


(4.30b) 


For a perfect gas they are identical to those derived by van Leer. The expression for the third 
component which satisfies continuity and symmetry conditions can be written most generally in the 
form 


•« = * 2ftTTIj £(« ± c ) 2 [ 2c ± (7 - i)**] 2 ± i^( u ± c ) ; 

± Pj-{ u ± c ) 2 ( u T c) 2 , 

4 c 


1(7-1) (7 — 1). 


(4.30c) 


where /3 is an arbitrary parameter. We again have a one-parameter family of flux splittings. 

For a perfect gas, the second term in vanishes. Van Leer’s condition of a vanishing eigenvalue 
for requires /? to vanish also, so that Ff^ t reduces to one term. One can readily show that the 
remaining eigenvalues are both of the proper sign for 1 < 7 < 3. For a real gas, for which 7^7, 
one requires at least two terms for F^ e . (The one-term solution of Montagne [83] actually violates 
the condition F = F + -f F~ .) Since 7 and 7 are both variable, it is imposible to obtain the vanishing 
eigenvalue condition identically throughout the velocity range for any choice of /3. 

One can demonstrate readily that for a thermally perfect gas, the two-term solution has one eigen- 
value that is of the wrong sign for the whole subsonic velocity range. Fortunately, the magnitude of 
the offending eigenvalue is extremely small, so that its effect on a numerical scheme is not noticeable. 
In view of this fact, and the general nature of the variation of 7 and 7, it is simplest to set /3 = 0 for a 
real gas. Further numerical experiments over a wide range of conditions are required to validate this 
approximation . 


4.3. Extension to Systems via the Local-Characteristic Approach 

In this section, the method sometimes referred to as the “local characteristic approach” (a gen- 
eralization of Roe’s approximate Riemann solver) in conjunction with TVD schemes is discussed. 
Consider a system of hyperbolic conservation laws of the form (4.1a) where U and F are vectors of 
m components. The idea of the local-characteristic approach is to extend the scalar TVD method 
to systems so that the resulting scheme is TVD for the “locally” frozen constant-coefficient system. 
The procedure is to define at each point a local system of characteristic variables W and to obtain a 
system of uncoupled scalar equations 

W t + KW X = 0, W^R~ l U (4.31a) 

A = diagfa*). (4.31b) 

The matrix R can be (4.6) for the Euler equations, and is a transformation matrix such that A is 
diagonal. One then applies the nonlinear scalar scheme to each of the scalar characteristic equations. 
The final form after transforming back to the original variables looks like the scalar case except there 
is coupling between the characteristic variables through the eigenvectors R. The numerical flux is of 
the form 

F h , = \{F i + F 1+l + fi, + 1 *, + il- (4.32) 

The matrix Rj+i is R evaluated at some symmetric average of Uj and Uj +1 . For example, Rj + i = 
R{ U3+ 2± ~ " ~ ) • Other approximate ways of obtaining a symmetric average are the Roe average as 
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discussed in sections 4.1 and 4.2, or those of Huang [30] and Carofano [81]. The /th element of 3 > j+ a 
denoted by <p l j + ± fo r the /th characteristic, has the same form as the scalar case except that Qj+i is 
replaced by a*. + i and A J + i is replaced by ai*. + i . Here are elements of <*j+i = R~*±(Uj+i -Uj). 
Specific forms of the <f> l J+ i f° r explicit and implicit methods will be discussed in sections 4.4 and 4.5 

The local-characteristic approach is more efficient than the exact [28] or approximate Riemann 
solvers of Osher and Solomon [76], and it provides a natural way to linearize the implicit TVD schemes 
[10,43]. The advantages of this approach as opposed to Davis’s simplified approach [8] or the Boris 
and Book approach [38] to systems are that (a) the current approach in effect uses scalar schemes on 
each characteristic field so that the limiter used need not be the same for each field (e.g., one can use 
a more compressive limiter for the linear fields and use a less compressive limiter for the nonlinear 
fields as in the numerical examples of references [10,44,84-86]); and (b) one can even use different 
schemes for different fields. For the one-dimensional Euler equations, the characteristic fields consist 
of two nonlinear fields u±.c and a linear field u. The contact discontinuities are associated with the 
linear fields. It has been shown [10,44,84-86] that the two different fields required different amounts of 
numerical dissipation (i.e., different limiters). Often the limiters that are designed for the linear field 
might give spurious oscillation or nonphysical solutions for the nonlinear field. Numerical examples 
concerning this aspect will be discussed in subsections 4.4.3 and 4.4.4, and in section 5. 


4.4. Description of the Explicit Numerical Algorithms and Examples 

The second-order in time and second- or third-order in space explicit-difference schemes considered 
here for both the MUSCL and the non-MUSCL approaches for the system case can be written in the 
same form as equation (4.9). 


4.4.1. The Non-MUSCL Approach 

The numerical flux functions for a non-MUSCL- type approach for both the upwind and 

symmetric TVD schemes [44,10] using the local-characteristic approach are given by equation (4.32). 

g 

Second-Order Symmetric TVD Scheme : The elements of the vector 3>y + i, denoted by (<^. + x ) , for a 
general second-order symmetric TVD scheme are 


= -X(a l . + ,) 2 Q 1 ^ 


3 + 


3+ o 


v+i- 




3 + h 


3 + *. 


(4.33a) 


The value a 1 . k is the characteristic speed a 1 evaluated at some symmetric average of Uj and Uj+\. 

3 + 2 

The function ip is an entropy correction to |z|. It can be the same as (3.18), which is repeated here 
for completeness: 




| z\ \z\ > Si 

(z 2 + 6 2 )/26! |*| < 6, • 


(4.33b) 


For the test problems containing only unsteady shocks to be shown, S\ is set to zero in most of 

the computations. Note that entropy-violating phenomena occur only for steady or nearly steady 

shocks. For steady-state problems containing strong shock waves, a proper control of the size of Si is 

very important, especially for hypersonic blunt-body flows. A discussion is given in section 5.7. The 

‘limiter’ function Q l , i, expressed in terms of characteristic variables, can be of the form 
3+2 
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g'. + 1 = minmod(a' _ i , a‘ i+ i ) + minmod(a'. + 1 , <*' + f ) - a' + 1 , 
Q‘ + i = minmod(a‘ a' a' ), 


— minmod 


^ a 3- k » i ’ ^ a i+-£ > o ( a i- k + a i+ # ) 


(4.33c) 

(4.33d) 

(4.33e) 


j+ £ ’ J+'f ’ 2 - 2 ’ 

The minmod function has the same meaning as in the scalar case. 

Second-Order Upwind TVD Scheme : The elements of the vector $ - + 1 for a second-order upwind TVD 

[7 ^ 

scheme , denoted by (<^ + i ) , originally developed by Harten and later modified and generalized by 
Yee [43,44,87], are 


W' + ±f = ff(a' + i)( 3 ' +1 + s') - <p(a‘ j+i + 7' + i)a‘ + i 


(4.34a) 


The function a(z ) — — Xz 2 \ and 




a »i = 0 ' 


(4.34b) 


Examples of the ‘limiter’ function g * can be expressed as 
9j = minmod (</_ ± , </ + x ) , 

o' = (“i+ * - * + l°j+ *“!■-* |) / (“S+ *+“$-*)• 

= {“i-i [K+i ) 2 + + “i+4 IK-i) 2 + **] 


( a y + l) + ( a ,_i) +26 2 


(4.34c) 

(4.34d) 

, (4.34e) 


9j = S •max|o,min(2|^. + i|,5 ■ a* _p,min(|c/ + i |, 25 • «' _i)|; 5 = sgn(aj. + p. (4.34f) 

The parameter 5 2 = 10 -7 in all of the calculations presented in sections 4.4.3 and 4.4.4 For a 1 + 
. . J+2 

<r ._ A = 0, g\ is set to zero in (4.34d). 

3 2 3 

4.4.2. The MUSCL Approach 

MUSCL Approach Using an Approximate Riemann Solver: The numerical flux function E . , i for a 

. • — — — — 3~r 2 

MUSCL type approach of an upwind scheme as described in Yee [22] using the local-characteristic 
approach can be expressed as 


'»■+* = + **+***+*]. 


(4.35a) 


where the elements of <&„• , i are 

J-r 2 


?‘ + i = 




(4.35b) 
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(4.35c) 


a 


J+: 




3+' 


3 + ' 


Here a 1 x are the eigenvalues and R + 1 is the eigenvector matrix of A, evaluated using a symmetric 

3~r 2 J 2 

average between Uf k and U*f k \ i.e., 

J ' 2 « 2 


a' +i =«‘(i/* 4 ,vf +i ), 


*** =*(ift*,itf + *)• 


J+ 2 
r£ 

’ J + 1 


(4.35d) 

(4.35e) 


However, there are options in applying the limiters for system cases. Namely, one can impose limiters 
on the conservative, primitive, or characteristic variables. For a first-order time-discretization, the 
simplest case is to define each of the elements of U R ’ L as in the scalar case (equation (3.38)). For a 
second-order time-discretization, in addition to the various options in imposing the limiter, another 
step is needed. For the moment, let us assume that the variables for imposing the limiter are W = 
(p,u,pe). Denote P and P~ 1 as the matrices such that U — PW and W — P~ l U . The vectors U R + L 

and U L i for a second-order in time, second-order in space MUSCL approach can be 

2 


Uf +i =P(W" +i + l -g,), (4.36a) 

V* +i =P(W% 1 i - ijy+,). (4.36b) 

Here gj is defined as in equations (4.34c)-(4.34f), except the arguments will be elements of (W" +l — W?) 
and (Wj - with = P - 1 (£/" + = ) where 

u" + ‘ = V? - \ ( F [ P(W " + |s,-)] - P [P(W? - \g,)\ | . (4.36c) 

One can define W = U the conservative variables or W — W the characteristic variables. A second- 
order in time but third-order in space scheme can be obtained by defining a different gj function (see 
equation (3.38) for a formula). 

MUSCL Approach Using Flux-Vector Splittings: The numerical flux *,•+* for either flux-vector split- 
tings, can be expressed as 


F ] + . = F+(V$ H ) + F-(Uf H ), (4.37) 

where F±(U L ! R ) are evaluated using either the generalized Steger and Warming splitting or the 

2 

generalized van Leer splitting. The vectors U R k and U L A are the same as in equation (4.36). 

° 3 + 2 j ‘ 2 


4.4.3. Comparative Study of TVD Schemes for Real gases [88] 

A detailed description of this study can be found in Montagne et al. [88]. Six one-dimensional 
shock-tube problems were considered. The left- and right-hand-side states of the initial conditions 
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for all six cases are tabulated in table 1. The cases have been ordered in the direction of increasing 
maximum Mach numbers encountered in the flow. 


Case 

State 

Density 

Pressure 

Temp. 

Energy 

Velocity 

Mach 

Case A 

Left 

0.0660 

9.84 10 4 

4390 

7.22 10 6 

0 

0.0 


Right 

0.0300 

1.50 10 4 

1378 

1.44 10 6 

0 

0.0 

Case B 

Left 

1.4000 

9.88 10 5 

2438 

2.22 10 6 

0 

0.0 


Right 

0.1400 

9.93 10 3 

2452 

2.24 10 6 

0 

0.0 

Case C 

Left 

1.2900 

1.00 10 5 

2718 

1.95 10 5 

0 

0.0 


Right 

0.0129 

1.00 10 4 

2627 

2.75 10 6 

0 

0.0 

Case D 

Left 

1.0000 

6.50 10 5 

2242 

2.00 10 6 

0 

0.0 


Right 

0.0100 

1.00 10 3 

346 

2.50 10 5 

0 

0.0 

Case E 

Left 

0.0100 

5.73 10 2 

199 

1.44 10 5 

2200 

7.8 


Right 

0.1400 

2.23 10 4 

546 

4.00 10 5 

0 

0.0 

Case F 

Left 

0.0100 

5.73 10 2 

199 

1.44 10 6 

4100 

14.5 


Right 

0.0100 

5.73 10 2 

199 

1.44 10 5 

-4100 

-14.5 


Table 1. Initial conditions for the test cases. 

The densities are given in kg/m 3 , the pressures in Newtonsjm 2 , the temperatures in 0 Kelvin, the 
internal energies in kg(m/ sec) 2 and the velocities in m/sec. 

The thermodynamic properties of equilibrium air were obtained from the curve fits of Srinivasan 
et al. [89]. These curve fits give analytic expressions for 7 in several ranges of density and inter- 
nal energy. The values of 7, and /c are then calculated from the derivatives of these analytical 
expressions. The numerical solutions were compared with an “exact solution” computed by solving 
the Rankine-Hugoniot jump conditions and integrating numerically the characteristic equations in 
the expansion fan. Figure 4.1 shows the distribution of Mach number and of the two quantities 7 
and 7 defined in (4.8). Not only the changes in their values but also the differences between them 
are indications of departure from the perfect-gas case. These differences do not necessarily occur at 
very high temperatures, but at intermediate temperatures when the vibration is excited and when the 
dissociation reactions start. Results for cases B, C and E are discussed here. 

The combination of Riemann solvers and of differencing algorithms considered above yields five dif- 
ferent schemes: a symmetric non-MUSCL scheme, an upwind non-MUSCL scheme, and three MUSCL- 
type schemes, depending on the Riemann solvers. The study provided a check on the validity of the 
extended formulas, since theoretical prediction of their properties appears to be difficult due to the 
non-analytic form of the state equation. Comparisons were done on the accuracy and the robustness 
of the methods. The six test cases chosen were intended to highlight the effect of the high ratios 
in pressure or density related to shocks, and the effect of departure from a perfect gas in the state 
equation. 

The five second-order explicit schemes tested are (a) the symmetric TVD scheme (4.33), (b) the 
upwind TVD scheme (4.34), (c) the upwind TVD scheme (4.35), (d) the generalized van Leer flux 
splitting (4.37) and (e) the generalized Steger and Warming flux splitting (4.37). Schemes (a) and 
(b) follow the non-MUSCL approach, while schemes (c)-(e) follow the MUSCL approach. The same 
approximate Riemann solver (local-characteristic approach, section 4.3) is used in the three schemes 
(a), (b) and (c). The limiter function used for each scheme remains the same for all the cases. Limiter 
(4.33e) is used for the symmetric TVD scheme (4.33), and limiter (4.34e) is used for the rest. The 
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time-step limit is expressed in terms of a CFL number related to the eigenvalues of the numerical 
fluxes. The CFL number is fixed at 0.9 in cases B and E. In case C, it has been fixed at 0.5 for 
the upwind non-MUSCL and MUSCL scheme and at 0.9 for the symmetric scheme. The actual CFL 
used for the flux splitting approaches is approximately 80% of the fixed CFL (see reference [88] for 
details). The number of discretization points is 141 in cases B and C, and only 81 for case E because 
the expansion fan is replaced by a shock. The time for stopping the computation has been chosen for 
each case in order to use the full computational domain. For a Ax normalized to 0.1 meter, these 
times, given in 10 -3 seconds, are t = 3.0 for case B, t — 5.0 for case C and t — 3.2 for case E. The 
comparison study can be divided into two aspects, one on the differencing algorithms, and the other 
on the Riemann solvers. Figures (4.2)-(4.4) show the perfect-gas computations and figures (4.5)-(4.7) 
show the real-gas computations. In all of the computations 5^ (in equation (3.18)) was set to zero for 
the non-flux-split approach. The vector W in (4.36) is set to (p,u,pe). Since the resolutions of shock 
and contact discontinuities do not always behave the same for the different variables, the computed 
solutions for 7,7, Mach number, energy, velocity and density are shown here for completeness. 

Comparison of the Differencing Techniques : Parts (a)-(c) of figures (4.2)-(4.7) provide a comparison 
of the symmetric TVD and the upwind TVD non-MUSCL scheme, and the MUSCL scheme with the 
same approximate Riemann solver for both perfect and real gases. 

The three techniques give almost the same results in general and the differences are similar to 
those found for a perfect gas. The greatest difference occurs in test case C. But this case happens 
to be already a difficult one when the same initial conditions on density and energy are applied to a 
perfect gas. The major differences are between the symmetric scheme and the two upwind schemes. 
Although the symmetric scheme is generally more diffusive at the contact discontinuities, the situation 
is reversed in case E where the main shock is almost stationary and the flow behind it has a very low 
velocity. Furthermore, this symmetric scheme yields more stable results in that case. The influence 
of the limiters is the same as for a perfect gas as summarized in references [10,44] and in sections 5 
and 6. A comparative study of flux limiters for case B will be discussed in the next subsection (4.4.4). 
This point, like some aspects related to the computational efficiency, is discussed more fully in [88]. 

The main difference in computational effort lies in the MUSCL and non-MUSCL approaches. The 
operations count between the non-MUSCL and MUSCL is within 30% for a perfect gas. However, 
due to extra evaluation in the curve fitting between the left and right states in an equilibrium real 
gas for the MUSCL formulation, additional computation is required for the MUSCL approach. The 
slight advantage of MUSCL over non-MUSCL is that MUSCL can be spatially third-order accurate. 
One word of caution is that experiences with the third-order case (( rj = 1/3 for equation (3.38))) 
do not show a very visible improvement over the second-order case. Part of the reason is that all 
TVD schemes reduce to first-order at points of extrema regardless of the order of accuracy in smooth 
regions. 

Comparison of the Riemann Solvers : Parts (c)-(e) of figures (4.2)-(4.7) compare the three Riemann 
solvers for the MUSCL scheme for both perfect and real gases. 

In the test conditions used, a comparison between the two classes of flux splittings showed little 
difference when equilibrium air was used. The generalized van Leer splitting yields a sharper capture 
of the shocks than the generalized Steger and Warming splitting. The resolution and operations count 
of the different approximate Riemann solvers [30], [81], [57] are very similar. The results obtained 
with the approximate Riemann solver are very similar to the ones obtained with the generalized 
van Leer splitting. Actually, the generalized van Leer splitting seems to be less sensitive to the state 
equation for the shock resolution while the approximate Riemann solver is more accurate at the contact 
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discontinuities. 

It is important to note that flux-vector splittings make use of the sound speed only, whereas ap- 
proximate Riemann solvers of the Roe-type make use of the thermodynamic dervatives \ and K °f 
equation (4.2). These thermodynamic derivatives put more stringent requirements on the curve fit 
that represents the thermodynamic properties of the gas. In this regard, the curve fits of Srinivasan 
et al. may be deficient for the approximate Riemann solver as can be seen from figure (4.1), case D. 
One probably needs more improved curve fits than those of reference [89] before a definite conclusion 
can be drawn about the accuracy of the different Riemann solvers and schemes. 

In conclusion, for the purpose of calculations in gas dynamics with equilibrium real gases, these 
numerical tests show that the simple extensions to a real gas for the flux-vector splitting or the 
approximate Riemann solver presented in this paper are valid. The main effect of using a real-gas 
equation of state is to exacerbate the problems of the methods for large discontinuities. Test case 
C is an example of such a situation. Similarity, it seems difficult to give a ranking of the methods. 
Depending on the case, each one presents some drawbacks or some advantages. The present results 
also indicate that the state equation does not have a very large effect on the general behavior of these 
methods for a wide range of flow conditions. 

None of the differences observed for the explicit versions seems to be decisive for the one-dimensional 
tests, but factors such as stability and computational efficiency need further investigation in multi- 
dimensional tests. The main differences between the methods lie in their versatility in extending to 
implicit methods with efficient solution procedures, especially for multidimensional steady-state com- 
putations. Preliminary study shows certain advantages of the approximate Riemann solver over the 
flux-vector splitting approaches (see section 4.5 for a discussion). 


4.4.4. Comparative Study of Flux Limiters for a 1-D Shock- Tube Problem 

The same five different explicit schemes discussed in the last subsection (4.4.3) were used to study 
the effect of different limiters on the accuracy of the schemes for both perfect and real gases. Because 
of the anticipated decrease in accuracy for real gases of the curve fits of Srinivasan et al. in providing 
the thermodynamic derivatives, only the result of case B is summarized and shown here. 

Figure (4.8) provides a perfect-gas comparison of the five schemes in conjunction with the different 
limiters (4.33) and (4.34) corresponding to the appropriate schemes. Figure (4.9) provides the same 
comparison for a real gas. A CFL = 0.9 was used for all computations. Parts (a) - (c) of figures (4.8) 
and (4.9) show the accuracy of Roe’s first-order scheme, Harten’s original modified-flux scheme (3.45) 
without the added artificial compression [6], and the author’s modification to (3.45) (i.e. (4.34a)- 
(4.34c)) respectively. One can see the dramatic improvement of scheme (4.34a)-(4.34c) over (3.45). 
Harten’s scheme (3.45) produced accuracy similar to that of Roe’s first-order method. 

Parts (d) - (m) of figures (4.8) and (4.9) summarize the non-MUSCL and MUSCL approaches for 
both perfect and real gases for all the second-order methods in section 4.4. In the case of a non- 
MUSCL approach, limiter (4.33e) is the most accurate among (4.33c) - (4.33e) for the symmetric 
TVD scheme (4.33). As for the upwind schemes, limiters (4.34d) and (4.34e) are very similar, whereas 
limiter (4.34f) gives very accurate results for contact discontinuities but is sometimes too compressive, 
thus causing slight oscillations in smooth regions for high Mach number cases (result not shown). A 
combination of limiters such as (4.34d) or 4.34e) for the nonlinear fields and (4.34f) for the linear 
field seems to be a good compromise. In the case of the MUSCL approach, only limiters (4.34c) and 
(4.34e) were studied. Between the two limiters, (4.34e) produces higher shock-resolution than (4.34c) 
(comparison not shown). 
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4.5. Description of the Implicit Numerical Algorithms and Examples 

The corresponding implicit scheme of equation (3.21) for the system case via the local-characteristic 
approach can be written as 


u ? +1 





= U? - (1 - 6) A 




(4.38) 


Here $ has the same meaning as in equation (3.21). The spatial accuracy of the scheme depends on 
the form of the numerical flux functions. Implicit methods via the flux-vector splitting approaches in 
conjunction with a MUSCL formulation can be found in references [61-63] and will not discussed here. 

Non-MUSCL Approach : The numerical flux function for a non-MUSCL type approach for both 

the upwind and symmetric TVD schemes has the same form as equation (4.13). The elements of the 

£ 

vector denoted by (^ +i ) for a general second-order symmetric TVD scheme are 


w$+*) =-^K+*)K+*-<?J + *). (4-39) 

where a 1 . A are elements of R~}i{Uj. fi — Uj). The function ip is the same as in the corresponding 

3 ' 2 3 ' 2 

scalar case. The limiter function Q l . x can be the same as (4.33). 




The elements of the vector + i denoted by (^. + i ) for a second-order implicit upwind TVD 
scheme, originally developed by Harten [6] and later modified and generalized by Yee [43,44,87], are 


. 1 / 1 


(<+*) = 2^K-+i)foJ+i + + ii+iK+4- 

The 7 1 . i and the limiter function g\ can be the same as (4.34). 

3 « o ^ 


(4.40) 


MUSCL Approach: The numerical flux function for a MUSCL-type approach of the implicit 

TVD upwind scheme has the same form as the explicit case except that the form of the U R ’^ is 

7+2 

slightly different. The values U R A and U L L are the same as (4.36) without the additional step 
(4.36c), i.e., setting 2 = Uj. 

A Conservative Linearized Form for Steady-State Applications : A conservative linearized form of 

(4.38) can be written as 


with 


i + o\{h; h 




E = t/ n+1 - U n , 


where 





n 


(4.41a) 


(4.41b) 

(4.41c) 
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Note that the matrix H x x in (4.41c) is different from the total enthalpy H in section 4.1. The 

2 

nonstandard notation n 

*?+** = \ (a h 1 E h1 - nj +1 s) (4.41d) 

is used, and for a first-order implicit operator, Q x k E can be 

2 

= fl J+ ,diag[0( a ' + i )]«;;,(£,+, - Ei). (4.41e) 

Here A,+ 1 , is the Jacobian of F evaluated at (j-j-l) and Ej = U* +1 -U”. For steady-state applications, 
one can simplify (4.41e) as 

n*. H E=M x I(E i+l -E j ). (4.41f) 

The scalar value .M x is 

M x - max i/j(a l . i), (4.41g) 

l J ' 2 

and / is the identity matrix. Note that (4. 4 If) involves scalar multiplications only. The solution using 
(4.41) is still second-order (or third-order) accurate after it reaches steady-state. A second-order form 
of (4.41e) for the non-MUSCL formulation can be found in references [10,43]. The nonconservative 
linearized form as discussed in the scalar case can be obtained similarly. See references [10,43] for 
details. 

The author would like to point out an important distinction between the flux-vector splittings 
and the local-characteristic approach for implicit methods. Unlike flux-vector splitting approaches, 
implicit methods employing the local-characteristic approach (non-MUSCL or MUSCL) with first- 
order implicit operators such as (4.41f)) do not require the Jacobian of the F ± fluxes. In many 
instances, the Jacobian of F± is relatively difficult to obtain. A similar difficulty applies to the MUSCL 
formulations via the local-characteristic approach if a second-order implicit operator is desired. 

To show the stability and accuracy of the implicit method for steady-state application, results for a 
quasi-one-dimensional divergent nozzle problem [90] are presented. Figure (4.10b) shows the converged 
density distribution after 25 steps at a CFL number of 10 6 using 20 equal grid spacings. The value 
in (3.18) is set to 0.125. The solid line is the exact solution and the diamonds are the computed 
solution. Only 14 points are plotted. The six points not shown on both ends of the z-axis are equal to 
the exact solution. The solution looks very much like the explicit TVD scheme (figure 4.10a) except 
it has a tremendous gain in efficiency. Figures (4.10c) and (4.10d) show the same computation with a 
classical shock-capturing method and the first-order Steger- Warming flux-vector splitting methods. 
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V. Extensions of Nonlinear Scalar TVD Schemes to Higher-Dimensional Nonlinear Sys- 
tems 

At the present stage of development, truly two-dimensional schemes are still in the research stage. 
The available theory is too complicated for practical applications. Here the usual approach of applying 
the one-dimensional scheme for each direction in multidimensional problems is adopted. Therefore, 
highly skewed nonorthogonal grid distributions should be avoided. As will be illustrated in later sub- 
sections, this methods of extending one-dimensional methods to multidimensional nonlinear systems 
is quite satisfactory. 

In order to preserve the original second-order time-accuracy, the extension of the nonlinear scalar 
second-order explicit schemes to multidimensional problems can be accomplished by a predictor- 
corrector, linear multistage or Strang-type of fractional-step [91] (time-splitting) method. All the 
time-accurate numerical examples illustrated later utilize the time-splitting and predictor-corrector 
methods. Therefore, only these two methods will be described here. For the implicit methods, only a 
simplified form will be briefly discussed, since extensions to other forms follow the same idea. Nonit- 
erative relaxation implicit methods other than alternating direction implicit (ADI) methods will not 
be discussed. Extensive work in the area of relaxation methods in conjunction with van Leer flux- 
vector splitting for perfect gases has been performed by Thomas, Walters and van Leer; see references 
[62,63,92] for details. All of the examples in higher dimensions employ coordinate transformations, 
and the schemes will be written in transformed coordinates. The formulation described here is valid 
for both two- and three-dimensional systems of conservation laws. Only the two-dimensional case will 
be described. For three-dimensional formulation, one only has to add an extra dimension and the 
corresponding numerical flux. 

Consider a two-dimensional system of hyperbolic conservation laws 


(W dF(U) dG(U) 
dt ^ dx dy 


= 0 . 


(5.1) 


Here U , F(U ) and G(U) are column vectors of m components. 

A generalized coordinate transformation of the form f == £(z,y) and rj — r)(x,y ) which maintains 
the strong conservation-law form of equation (5.1) is given by 


dU dF(U) dG(U) 

at d( dr, 


(5.2) 


where U — U/ J , F — (£ X F + £ y G)/J , G — (r) x F + r) y G)/J , and J = £ x r) y - £ y r) x , the Jacobian 
transformation. Let A = dF/dU and B — dG/dU. Then the Jacobians A and B of F and G can be 
written as 

A = (£ X A + £ y B) (5.3a) 

B = (r} x A + r) y B). (5.3b) 

Let the eigenvalues of A be. (a|,a|, ...,a^) and the eigenvalues of B be (a*,a^, ...,a^). Denote R ^ 

and R n as the matrices whose columns are eigenvectors of A and B , and denote R^ 1 and i?" 1 as the 
inverses of R $ and R n . 

Let the grid spacing be denoted by A£ and A r) such that f = jA£ and rj = kArj. Let U j+ i k 
denote some symmetric average of U^ k and Uj+ ljk (for example, U j+ i k = 0.5 * (U 3+ itk + Uj.k), or 
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Roe’s average). Let a^. + 1 , Rj+i, R~+i denote the quantities a^, R R J 1 related to A evaluated at 
Similarly, let a l k+k , A denote the quantities a^, R^, R~ l related to B evaluated 

at 

In all of the two-dimensional calculations 


a 


2 


^-1 U 3 + l,k ~ Uj,k 

^ 2 0.5 * (*/j+i J fc H - 


is the difference of the characteristic variables in the local ^-direction, and 


(5.4a) 


oc k+ x = R 


-l 


U 


j.fc+i 


Ui 


fc +2 0.5 * (J jt fc+i + Jj )k ) 


(5.4b) 


is the difference of the characteristic variables in the local ^-direction. The symbol J 3tk is the Jacobian 
transformation evaluated at (j'A£, fcA*?). The averaged Jacobians are used here in order to preserve 
the free stream. 


5.1. Description of the Explicit Numerical Algorithms and Examples 

Fractional-Step Method: By using the Strang-type of fractional-step (time-splitting [91]) method, the 
one-dimensional second-order TVD scheme in section 4.4.1 can be implemented for the two-dimensional 
system (5.2) via the local-characteristic approach as follows: 


where 


frn- 1-2 
U 3,k 




= U U 


u ,V 

— (, 
A r) V 


h nh ph ph/2frn 

u j,k> 

(5.5a) 


(5.5b) 

T* P«* \ 

*3,k+\ G j,k-U 

(5.5c) 


with h — At. The functions Fj + i k and are the numerical fluxes in the f- and ^-directions 

evaluated at (j + |,fc) and (j, k + 2 ), respectively. Typically, Fj+i k can be expressed as 

Fj+±,k = 2 (Fj> k + Fj+i,k + Rj+iQj+i)- ( 55d ) 

Here Rj+x is the eigenvector matrix for dF/dU evaluated at some symmetric average of Uj ik and 

Uj + i t k- Similarly, one can define the numerical flux G k+ ± in this manner. The function in the 

^-direction for either the symmetric or upwind TVD schemes can be expressed in the same way as 
equations (4.33) and (4.34) with the appropriate eigenvalues and defined above. See reference 

[93,24] for details. 


Predictor-Corrector Method : To preserve the original second-order time- accuracy, besides using the 
time-splitting approach, one can employ a predictor-corrector method similar to the one discussed in 
the nonlinear scalar case. With the same notation, a formal extension of the scalar explicit second- 
order TVD method in predictor-corrector form (via the local-characteristic approach for the nonlinear 
hyperbolic system (5.2)) can be written as 
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(5.6a) 
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3 + l,fc 
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A T) 


G i* 


_ r (!) 
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(5.6b) 


^ +1 


y i,l’ + 


- r {2 \^ [2 \ 

3 2 3~\~ 2 3 2 3 2 
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F 



(2) 




(5.6c) 


Here the superscripts (1) and (2) designate the values of the function evaluated at the intermediate 
solutions and The /th element (f > 1 . A of the vector . i has the same form as the scalar 

3 ~ r 2 J ' 2 

case (3.59d) except the a j+ 1 is replaced by a^. +i and the Aj + i is replaced by of equation (5.4a). 

See sections 3.3, 3.6, 3.7 for the discussion of equation (5.6c) and the usage of the scheme (5.6) in 
general. 


5.2. Time- Accurate Computations by Explicit Methods 

Four time-accurate computations for a perfect gas are described here. All of the descriptions are 
summaries of the published or yet to be published works of authors indicated in the subsection head- 
ings. The first three computations are external flows and the last one is an internal flow problem. 
Experimental data were available for three of the problems. In all of these three problems, good 
agreement between the computations and the experimental data were observed. 


5.2.1. Shock Wave Diffraction From an Obstacle (Young and Yee [85]) 

There is a continuing interest in determining the diffraction loadings imposed on a stationary object 
during a blast-wave encounter, since this knowledge is important in designing the structure to survive 
such an event. The problem is especially suitable for numerical simulation, since experimental setups 
for such studies are usually very expensive and sometimes impractical. 

A generic configuration of a class of objects of interest is shown in figure (5.1). The configuration 
has a broad base to maximize stability. To reduce the drag force acting on the body, the top is rounded 
off to minimize vortex shedding and flow separation. The objective is to ensure that the downward 
force is much larger than the lateral force generated on the body during a blast-wave encounter so 
that the object would suffer only minimal lateral motion and would not tip over. 

Computational Domain and Grid Generation : For illustration purposes, a wedge angle of 40° and a 
rounded top with a radius of curvature of 0.17 times the base width was chosen for the current study. 
Different wedge angles and top curvatures were also computed but are not reported here. The flow 
features depend strongly on the wedge angle and top shape. The grids used were generated by a 
generalized Schwarz-Christoffel transformation for curved surfaces, and are orthogonal except at both 
ends of the body. Figure (5.2) shows the computational domain and grid distribution. A CFL number 
of 0.99 and a normal incident shock of Mach 2.0 in a perfect gas with 7 = 1.4 and £1 = 0 (t f>(z) of 
equation (3.18)) were used in the computations. Roe’s average was used to evaluate and a J+ x. 

Figure (5.3) shows a selected sequence of the diffraction process using the time-split second-order 
upwind scheme (4.34a) together with limiter (4.34d) for the nonlinear fields and (4.34f) for the linear 
fields. Limiter (4.34c) is the most dissipative limiter of the three limiters (4.34c)-(4.34f) and the results 
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are not shown here. Figure (5.3a) shows the distinct formation of the triple point, with the Mach stem 
and the contact surface emanating from the triple point. 

Figures (5.3b)-(5.3d) show the diffraction process as the shock wave traverses the rounded top. The 
triple point and the contact surface remain very distinct and the Mach stem has evolved into a curved 
diffracted shock. It remains almost perpendicular to the body surface at the impingement point. The 
diffracted shock is actually travelling slightly ahead of the incident shock. A careful examination of 
figures (5.3d)-(5.3f) reveals that the contact surface is beginning to roll up to form vortices. Figure 
(5.3g) shows the instant when the diffracted shock leaves the trailing edge of the body. Figure (5.3h) 
shows the fine resolution of the captured reflected shock at the trailing edge of the body and the 
emergence of a Mach stem located near x = 0.8. 

Figure (5.4) shows the same sequence of the diffraction process computed using the predictor- 
corrector symmetric TVD scheme (5.6) together with limiter (4.33e). The shock resolution is almost 
identical to that of figure (5.3) except with slightly more diffusive slipstreams. Figure (5.5) shows 
the density contours comparing the upwind scheme (4.34a) using limiters (4.34d,f) with the predictor- 
corrector symmetric scheme (5.6) using limiters (4.33c)-(4.33e). The cut off appearance of the incident 
shocks indicates that the incident shock has travelled slightly beyond the right edge of the computa- 
tional grid. Judging from the density contour plots of figures (5.4) and (5.5), it would appear that 
the result using (4.33e) is comparable to that of the upwind TVD scheme with limiter (4.34d) for the 
nonlinear fields and (4.34f) for the linear fields, although the flow field is slightly smoother than the 
upwind method away from discontinuities. Results from additional numerical experiments not shown 
here also indicate that the shock resolution of limiter (4.33e) is slightly better than the upwind scheme 
with limiter (4.34c). Considering that the predictor-corrector variant requires a smaller operation 
count and allows a larger time-step, it offers a very attractive alternative to the upwind TVD scheme. 

A comparison was also made between the time-splitting symmetric TVD scheme and the predictor- 
corrector scheme (5.6) using the same limiters ( (4.33c) - (4.33e) ). It was found that the predictor- 
corrector formulation produced slightly sharper shock resolution than the time-split symmetric form. 

Numerical Results by the Explicit MacCormack Method: Using a Courant number of 0.99 and the 
same fine grid of figure (5.2) for the explicit MacCormack method with a smoothing coefficient of 
0.2 resulted in negative pressure even before the incident shock reached the top of the body and even- 
tually led to program abortion. The Courant number was reduced to 0.5 and the negative pressure 
did not show up until the incident shock reached x = 0.182. A sequence of the results by the explicit 
MacCormack method is shown in figure (5.6). Figure (5.6d) shows the instant when the negative 
pressure first appeared. It is observed from the thickness of the shocks and oscillations associated 
with them that the shock-capturing resolution of the explicit MacCormack method is inferior to the 
TVD schemes. 


5.2.2. Shock Wave Diffraction from a Cylinder 

In Yee and Kutler [93], a time-split form of the modified-flux approach with added artificial compres- 
sion was used to simulate a planar moving shock wave impinging on a circular cylinder. The results of 
Yee and Kutler reveals the need of a better flux limiter for capturing contact discontinuities, and the 
need of fine enough grid resolution. Recently, V.Y.C. Young of Martin Marietta recomputed the same 
problem with an incident shock of Mach 2.81 on a finer grid using the predictor-corrector symmetric 
TVD scheme with limiter (4.33e). Roe’s average was used to evaluate R J + i and aq + ± , and 8i is 
set to zero in i/>( z ) (equation (3.18)). This particular Mach number was chosen because experimental 
Schlieren photographs were available from A. Bryson of Stanford University. The detailed description 
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of the shock-tube experiment can be found in Bryson and Gross [94] . 

Due to the symmetry of the problem, only the top half-symmetry plane consisting of the semi- 
annular region between the radii of 0.5 and 3.0 was considered. In the angular direction, 362 rays 
with a uniform angular spacing of half a degree were used. The first and the second rays straddled 
the symmetry plane, and similarly for the last and the next to last rays. In the radial direction, 202 
grid points were distributed with an exponential stretching. The first radial grid point was imbedded 
into the body as the image point, with the second radial grid point serving as the body grid point. 

A Mach 2.81 normal shock wave was initially located at x = —0.75. The computation is stopped 
at approximately the same position as the Schlieren photograph (x « 2). The results in figure (5.7) 
show that a fairly detailed flow structure was obtained, especially in the wake region. The locations 
of the contact discontinuities emanating from the triple points were clearly captured. The contact 
discontinuities close to the centerline inside the wake region were also captured. The formation of 
vortices could also be spotted inside the wake. Overall, the result closely duplicated the experimen- 
tal Schlieren photograph. Because the computation is inviscid, the locations of the discontinuities 
appeared to differ slightly from those of the experiment. 

5.2.3. Complex Shock Reflections Prom Airfoils at High Angle of Attack (Moon and Yee 
[84]) 

An interesting shock-tube experiment was conducted by Mandella [95] and Mandella and Bershader 
[96]. A schematic sketch of the experimental apparatus and the flow field are shown in figure (5.8). 
After the diaphram ruptures, a planar-shock wave of Mach number M s = 2 forms and travels down the 
shock tube. The shock tube has a 5cm x 5cm inner cross section and its end is open to the ambient 
atmosphere. In order to keep the shock wave and its induced flow two-dimensional, two parallel plexi- 
glass plates are attached to the end of shock tube. An NACA 0018 airfoil is mounted between the 
plates at an angle of attack of 30°. The airfoil is located at a distance of 3 1/3 times the shock tube 
height away from the shock tube end. The planar shock starts to diffract as soon as it leaves the 
shock tube, and forms a curved-shock wave which finally impinges on the airfoil. The curved-shock 
wave loses its energy in time as it diffracts and the Mach number decreases to approximately 1.5 at 
the instant of impingement. 

Due to the positive angle of attack of the airfoil, an interesting feature of shock reflection is observed. 
There is a transition from a short moment of compression to expansion on the upper surface, and 
compression on the lower surface. Along the upper surface, the shock reflection retains the regular 
reflection up to the compressive region and then makes a rapid transition from regular to Mach 
reflection. A single Mach reflection forms a triple point from which the Mach stem, contact surface, 
and reflected shock emanate. Also a vortex starts to form and grows very slowly at the upper nose due 
to the sudden but mild expansion. Meanwhile, the transition process develops gradually but rather 
strongly from a regular to a Mach reflection on the compressive lower surface. Eventually, the Mach 
stem developed on the lower side wraps around the trailing edge of the airfoil and a vortex is generated 
due to the sudden strong expansion and the sharp trailing edge. 

At the exit of the shock tube, another interesting and equally complex flow pattern is observed in 
the early stages of shock diffraction. The structure of flow fields behind the diffracting shock behaves 
in a pseudo-stationary way on a 90° convex corner until the Mach wave hits the plane of symmetry. 
A slipstream emanating from the convex corner rolls into the vortex-spiral. The reflected Mach wave, 
the incident planar shock, and its diffracted part form a triple point from which the contact surface 
emanates and ends at the slipstream encircling the vortex. The second shock appears and forms a 
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triangular configuration with the tail of the Prandtl-Meyer fan (the terminator) and the slipstream. 
After the Mach wave hits the plane of symmetry, most of the features in the flow field retain their 
structures for some moments. The second shock grows with time and forms into a Mach disk after 
hitting the plane of symmetry. 

In order to properly simulate the experiment, the curved-shock solution before it reaches the airfoil 
is required as an initial condition. Also time-dependent boundary conditions along the outer boundary 
are required. These requirements are met by simulating the whole experimental region from the shock 
tube to the airfoil. The computational domain would have to cover a large area with a reasonably 
fine grid throughout the entire domain in order to capture the complex shock structure as time 
evolved. To avoid unnecessary computations, to confine the physical domain to be smaller, and to 
make use of body-fitted coordinates around the airfoil, the numerical simulation is broken into two 
parts: (a) simulation of the shock tube that contains a 90° convex corner (with a computational 
domain extended to cover the region where the airfoil is supposed to be located but without the 
airfoil), and (b) simulation of a time-dependent curved shock interacting with an airfoil. Figures (5.9) 
and (5.10) show the computational domain for the two parts. The flow field solutions obtained in part 
(a) are stored and used with bilinear interpolation to obtain the necessary initial and time-dependent 
boundary conditions for the airfoil simulation. The time-dependent boundary conditions were obtained 
assuming that the presence of the airfoil had no effect on the flow at the outer boundaries. 

The numerical simulation of the experiment is also conducted in a simpler way. The curved moving 
shock impinging on the airfoil is modelled as a constant-velocity, planar incident shock with approxi- 
mately the same incident Mach number. Since the shock is planar, the implementation of the boundary 
condition along the outer boundary is easily handled with the known solutions. The numerical simu- 
lations of the curved- and planar-shock waves interacting with the airfoil will be discussed shortly. A 
better understanding of the difference in shock structure for the planar- and curved-shock simulations 
is useful for the vortex-blade interaction parametric study where a planar-shock wave is preferred. 
The second-order explicit upwind and symmetric TVD schemes with the Strang-type of time-splitting 
were applied to these cbmplex physical problems solving the Euler equations of gas dynamics. Roe’s 
average was used to evaluate Rj+i and a J+ 1 , and was set equal to 0 in all computations. 

Shock Diffraction on a 90° Convex Corner : Figure (5.11) shows three sequential instants of a diffract- 
ing shock wave of Mach 2 on a 90° convex corner. The time intervals are 50 //sec and 100 //sec. 
Computed density contours from the solution obtained with the second-order explicit upwind TVD 
scheme (4.34a) together with limiter (4.34f) for the linear fields and (4.34c) for the nonlinear fields 
at a CFL number of 0.5 are shown in the right-hand column. The computed results agree quite well 
with the interferograms at each corresponding instant. The interferograms shown in the left column 
are missing in each corner of the photographs because of the configuration of the test-section window. 
The diffracting shock and the second shock are well captured within three to four grid points. The 
thickness of the shocks in the region for x < 2.2 is twice that of those in the rest of the region because 
the grid spacing is 2 times coarser in that area. The contact surface, slipstream, and the Prandtl- 
Meyer fan are also well captured. Their interaction with the vortex is shown in the first two instants. 
In the current case of a Mach 2 shock wave, the vortex core is not clearly defined. In fact, as the 
Mach number of the shock wave increases, the vortex core becomes more diffuse [97]. This appears in 
the numerical computations as well as in the interferograms. Numerical experiments carried out for 
this case show that a CFL number greater than 0.5 resulted in instability at the vortex core. Further 
study on grid resolution and the proper use of limiters for vortex- type flow fields is needed. The third 
instant (in which the test-section window is shifted downstream) shows a Mach disk formed after the 
second shock hits the plane of symmetry. The computed result also shows the well-captured Mach 
disk within 3-4 grid points following the contact surface. 
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Curved Shock Interaction with an Airfoil : Figure (5.12) shows the numerical results at six sequential 
instants while the traveling curved shock is progressing over a NACA 0018 airfoil at an angle of attack 
a — 30°. The time intervals are 10 /i sec except for the second interval which is 20 //sec. Computed 
density contours are shown for the solution obtained using the second-order explicit upwind TVD 
scheme (4.34a) together with limiter (4.34f) for the linear fields and (4.34d) for the nonlinear fields 
with a CFL number of 0.98. It should be noted that if the limiter given in equation (4.34f) is used for 
all of the fields, the solution diverges at a very early stage of the interaction. The authors conjectured 
that limiter (4.34f) is too compressive for the nonlinear fields. It was found that in general, the use of 
(4.34d) for the nonlinear fields produced sharper shock resolution than (4.34c). The use of (4.34f) for 
the linear fields produced sharper contact discontiniuties than (4.34d). However, the use of (4.34d) for 
the nonlinear fields for the computation of shock diffraction on a 90° convex corner will give divergent 
solutions near the corner and the vortex core regions. 

The Mach number of the time-dependent curved-shock wave is approximately 1.5 at the moment of 
impingement. A strong vortex and the dense contours of density around it appear in figure (5.12e), 
since the Mach stem turns over an approximately 180° convex corner. The most inclusive features of 
the shock interaction with the airfoil are contained in figure (5.12f). At this instant, the foot of the 
Mach stem proceeds upstream and collides with the incident moving shock on the upper surface. The 
experimental interferogram and computed density and pressure contour plots at this time are shown in 
figure (5.13). It should be pointed out that the two little bumps on both the upper and lower surfaces 
two thirds of the way along the airfoil and the very thin layer wrapping all around the airfoil observed 
in the interferogram result from vibration of the experimental setup. Also, because the airfoil used in 
the experiment is a hand-made NACA 0018 airfoil, a slight difference exists between the experimental 
and computed configurations. The incident and reflected shocks, Mach stems, and contact surface on 
the lower surface are captured within three grid points. Also the vortices at the upper nose and the 
trailing edge of the airfoil are well captured. On the upper surface, the experimental results show a 
weak contact surface, which is not captured in the simulation. By increasing the grid resolution around 
that region, this slip surface is also expected to be captured. As expected, the contact surfaces cannot 
be seen in the pressure contour plot because the pressure is continuous across the contact surfaces. 


Planar Shock Interaction with an Airfoil : For comparison, a planar incident shock of Mach 1.5 at the 
moment of impingement was used to approximately model the curved shock to study the difference 
in the computed flow field structures. Figure (5.14) shows a schematic representation of the physical 
plane with its boundaries and initial conditions. The same C-grid as in the previous case is used. The 
initial condition is implemented in the same manner as in the shock-tube section. Along the outer 
boundary, the analytic boundary condition is used to trace the location of the moving planar shock 
as a function of time. The other boundary conditions (along the body and the wake region) are also 
implemented in the same manner as was used for numerical simulation of the curved incident shock. 
Figure (5.15) shows four sequential instants of the interferograms compared with the density contour 
plots of the numerical results obtained for both the curved and planar shocks. The time interval 
between the first and second interferograms is 20 //sec and the rest are 10 //sec apart. For both cases 
the numerical results are computed using the second-order explicit upwind TVD scheme at a CFL 
number of 0.98. The curved-shock solutions agree quite well with the interferograms. The planar- 
shock solutions also compare favorably with the experiment in terms of the shape and location of the 
discontinuities. There is, however, a slight difference in density level between the two simulations. 
The same format of comparison is shown in figure (5.16) using the second-order explicit symmetric 
TVD scheme (4.33a,e) for both the curved and planar shocks. The shock resolution is similar to the 
upwind TVD scheme except that it is slightly more diffusive. The shock resolution for the symmetric 
TVD scheme (4.33a,e) is, however, sharper than (4.33a,c) and (4.33a,d), and the upwind TVD scheme 


54 



(4.33a,c). 


The numerical experiments show that both upwind and symmetric TVD schemes are quite stable and 
accurate even for higher Mach number shocks. The study also shows that for higher- Mach- number 
cases the symmetric TVD scheme is less sensitive to the numerical boundary condition treatment 
than the upwind scheme. For M s > 10, one has to use a characteristic type of numerical boundary 
condition for the upwind TVD scheme. Aside from the difference in numerical boundary treatments, 
one advantage of symmetric TVD schemes over upwind TVD schemes is that the symmetric TVD 
schemes require less CPU time than the upwind schemes. Computations with the symmetric TVD 
scheme (4.33a,e) at a CFL number of 0.98 for a planar incident shock of M s — 20 and a — 30° is shown 
in figure (5.17). The grid used is the same as for the previous cases. It is interesting to note that in 
contrast to the previous case, the vortex formation due to the sudden expansion is not noticeable and 
the reflection shock is swept around the body. 

5.2.4. Shock Propagation in a channel with 90° Bends (T. Aki [86]) 

This computation was performed by T. Aki of the National Aerospace Laboratory, Tokyo, Japan. A 
detailed description of the physical problem and numerical computations are reported in his original 
paper [86], 

A shock wave transmitted through a bend in a channel modifies its shape as it travels around the 
bend, although its initial shape is planar and the bend curves smoothly. The modification of the shape 
results in an irregular force acting on the bent wall and nonuniform flow behind the transmitted shock. 

The physical processes taking place during the transmission were little known and only qualitative 
arguments based on the investigations of the shock processes for isolated concave and convex corners 
had been given. An experimental investigation of the shock process through 90° bends was performed 
by Takayama et al. in 1977. One can find the details of the experimental setup in reference [98]. 

A Mach reflection established on the outer wall communicates to a diffracted wave around the inner 
wall. The triple point moves successively toward the inner wall. The reflected wave emanating from 
the triple point interacts with a contact surface appearing behind the diffracted wave and eventually 
hits the inner wall and then proceeds along it. 

The photograph (figure (5.18)) is an infinite fringe interferogram at this stage. The incident shock 
Mach number (at the entrance of the bend) is 2.2 and the test gas is air. The inner radius of curvature 
of the testing bend is 80 mm and its width is 40 mm. Therefore, the nondimensional inner curvature 
of the bend based on the width is 2. Takayama et al. concluded that the curvature of the bend must 
be equal to or greater than 2 in order to obtain a recovery of the planar shock after transmission 
under their experimental conditions. Generally speaking, the higher the incident shock Mach number 
is, the larger the curvature of the bend needed for a stable shock transmission; i.e. a recovery of planar 
shock front. After the stage shown in the photograph (figure (5.18)), the Mach reflected and diffracted 
shocks merge into one wave and one can observe recovering in this case a planar shape within the 
bend itself. 

The numerical boundary condition treatments are as follows: Let j = 1 be the inflow boundary. 
If the flow behind the incident shock is supersonic, then the variables on j = 1 are fixed with those 
obtained from the moving shock relations. If it is subsonic, then a procedure similar to that on the 
wall is employed. The variables on j = 1 are updated by using the Riemann invariants for the inflow 
and outgoing characteristics. The Riemann invariants for the inflow characteristics are postulated as 
those located far upstream; i.e. those behind the incident shock. At the outflow, the computation 
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terminates when the transmitted shock front or part of it arrives at the outmost grid plane. 

At the walls, let k — 1 be the grid plane on the inner wall. Initially, the first-order upwind scheme is 
used to evaluate U j at the inner wall. Since this scheme contains only information for the left-running 
characteristics, Uf must be updated to t/” +1 taking into account the effect of the right-running or 
reflected characteristic. The Riemann invariant for this wave, the entropy equation, and the boundary 
condition (vanishing normal velocity) are sufficient for updating . Treatment on the outer wall at 
k — 121 is similar to reversing the characteristics to be considered. 

In the computation, a total grid size of 581 X 121 is used with 450 x 121 of the grid located in the 
bent section. In the case of the time-split upwind TVD schemes using limiter (4.34f), a 400 x 121 
grid in the bent section was used. A CFL number of 0.98, 8\ = 0.1 and Roe’s average were used 
in all computations. Figures (5.18) and (5.19) show the pressure and density contours compared 
with experimental data for both the time-split upwind and symmetric TVD schemes using three 
different limiters. For the upwind TVD schemes, limiters (4.34d), (4.34d,f) and (4.34f) were used, 
corresponding to figure (5.18a, b,c). For the symmetric TVD schemes, limiters (4.33c), (4.33d) and 
(4.33e) were used corresponding to figure (5.19a,b,c). The overall performance compares very favorably 
with experimental data. The symmetric TVD scheme compares closely with the upwind TVD scheme, 
and with better accuracy than in the airfoil problem discussed earlier [84]. If one sets — 0 (equation 
(4.33b)), a slightly better shock-resolution is expected (e.g. see reference [99,100]). The result using 
6\ — 0 and the result using the predictor-corrector version (5.6) are under investigation and will be 
reported in Aki’s final paper. 

The main difference in accuracy at the walls for the four time-accurate calculations here is that 
each of these computations used slightly different outflow and wall numerical boundary conditions. 
Moreover, a different computer, different computer implementation, and a different density in grid 
spacing were used in these examples. 


5.3 Description of the Implicit Numerical Algorithms and Examples 


A one-parameter family of explicit and implicit TVD schemes for the two-dimensional system (5.1) 
can be written as 


^ +1 + + - ® X -*1 


= ^ - *"-* ,J - - d l> 




(5.7) 


Here & 3 + x and have the same form as in the implicit method defined in section 4.5 

The solution procedures for steady-state calculations in two-dimensional and three-dimensional 
Euler and Navier-Stokes equations are as follows. For the explicit operator, the convection terms are 
discretized by TVD schemes, and the diffusion terms are approximated by a central-difference method. 
For the implicit operator, a linearized conservative delta form can be constructed. For efficiency, a 
spatially first-order implicit operator, as in equation (4.41), was employed in most of the experiments. 

For steady-state applications, the resulting linearized delta form can be solved by some appropri- 
ate relaxation method other than ADI. This is the direction of current research. However, only a 
conservative noniterative ADI form [43] for the homogeneous PDE’s will be described below. For 
steady-state applications, the numerical solution is independent of the time-step. The implicit opera- 
tor has a regular block-tridiagonal structure and the resulting block tridiagonal matrix is diagonally 
dominant. One can modify a standard central difference classical shock-capturing code by simply 
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changing the conventional numerical dissipation term into the one designed for the TVD scheme. The 
only difference in computation is that the current scheme requires a more elaborate dissipation term 
for the explicit operator. No extra computation is required for the implicit operator. 


A Conservative Linearized ADI Form for Steady-State Applications : A conservative linearized ADI 
form of equation (5.7) used mainly for steady-state applications as described in detail in reference 
[10,43], can be written as 
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Here are (5.3) evaluated at (j+l,k) and 1). All of the inviscid calculations shown 

in this paper use (5.8) for steady-state applications. Other linearized forms suitable for time-accurate 
calculations are reported in references [10,43]. 


5.4. Steady-State Computations by Implicit Methods 

The numerical experiments were mainly performed on a NACA 0012 airfoil using the local- 
characteristic approach with = 0.125. For subsonic to low supersonic flows, the resolution of 
the shock waves was found to be quite insensitive to 0.1 < Si < 0.125 and a constant value seems to be 
sufficient. However for hypersonic flows, especially for blunt-body flows, a constant Si was found to be 
insufficient, and a variable Si depending on the spectral radius of the Jacobian matrices of the fluxes is 
needed. Moreover, a proper choice of the entropy parameter S\ for higher Mach number flows not only 
helps in preventing nonphysical solutions but can act, in some sense, as a control in the convergence 
rate and in the sharpness of shocks and slip surfaces (or shear layer in viscous flows). The smaller the 
nonzero being used, the slower is the convergence rate. The larger the S x being used, the larger is 
the numerical dissipation being added. See section 5.7 for numerical examples. 
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Generally, for inviscid time-accurate calculations, upwind TVD schemes produce sharper shocks 
than symmetric TVD schemes [8]. For the current implicit symmetric TVD scheme with limiter 
(4.33c) or (4.33d), this seems to be not the case. The symmetric method appeared to produce almost 
identical results as those from an upwind TVD scheme (4.34a)-(4.34c). 

Numerical studies also show that there is no difference in resolution in using limiter (4.33c) or 
(4.33d) for the symmetric TVD scheme. Limiter (4.33e) produces slightly sharper shocks than (4.33c) 
and (4.33d). This conclusion was based on a numerical study for flow-field conditions ranging from 
subcritical to transonic and supersonic for the NACA 0012 airfoil. Also, since these test cases consist 
of shock waves only, the same limiter was used for all characteristic fields. Figures (5.20) and (5.21) 
show a comparison of the current method using limiter (4.33c) with the upwind scheme (5.1) for 
two inviscid steady-state airfoil calculations. The two solutions are almost indistinquishable. For the 
current calculations, the upwind TVD scheme requires approximately 35% more CPU time than the 
symmetric TVD scheme on the CrayX-MP computer. 

Figures (5.22) and (5.23) show an inviscid comparison of the symmetric TVD scheme with the widely 
distributed computer code ARC2D, version 150 [101] . The free stream Mach numbers are Moo = 1.2 
and 1.8, and the angle of attack is a = 7°. The pressure coefficient distributions (not shown) are 
identical between the two methods and yet the flow field appears very different. The symmetric TVD 
scheme gives a very well-ordered flow structure and can still capture the shocks with a coarse grid, 
especially near the trailing edge of the airfoil. On the other hand, the ARC2D code did rather poorly. 
The ARC2D, version 150 computer code is based on the Beam and Warming ADI algorithm [17], but 
uses a mixture of second and fourth-order numerical dissipation terms. These numerical dissipation 
terms contain adjustable parameters. The values of the parameters on figures (5.22) and (5.23) are 
the same value as suggested in reference [101]. Other values of the parameters besides the one used in 
reference [101] were also studied. What is shown here is representative of the performance of ARC2D 
for this range of Mach numbers and angles of attack. For subsonic and transonic flow regimes the main 
advantage of TVD schemes over ARC2D for steady-state calculations is that one can capture the shock 
in one to two grid points as oppose to three to four. The flow away from the shock looks very much 
like that calculated by ARC2D. Note that in general, for two-dimensional blast wave calculations, the 
symmetric schemes usually can capture the shock in two to three grid points, but the slip surfaces are 
slightly more diffusive than the upwind TVD schemes as discussed in section 5.1. 

The same problem was studied for the upwind TVD scheme, and the results and convergence rates 
were found to be almost identical to those for the symmetric TVD scheme. For figures (5.22) and 
(5.23), a residual of 10“ 12 can be reached at about 400-600 steps. ARC2D, however, required only 
200-300 steps to converge to the same residual. 

Steady-state computations for higher Mach number flows in the hypersonic regime are currently 
under investigation. Preliminary study shows that the ADI form of the implicit scheme is very sen- 
sitive to initial conditions and numerical boundary condition treatments for blunt-body flows. The 
convergence rate is lower than subsonic and supersonic flows. 


5.5. A Thin-Layer Navier-Stokes Calculation 

For steady-state application, a simple algorithm utilizing the TVD scheme for the Navier-Stokes 
equations is to difference the hyperbolic terms the same way as before, and then central difference 
the viscous terms. The final algorithm is the same as equation (5.8) except that the spatial central 
differencing of the viscous term is added to the right-hand side of equation (5.8). The numerical 
solution shown below illustrates that this algorithm produces a fairly good solution for the case of an 
RAE 2822 airfoil calculation. Figure (5.24) is an example of the viscous case for the RAE 2822 airfoil 
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using the implicit upwind TVD scheme using limiter (4.34c). The thin-layer Navier-Stokes equations 
with the algebraic turbulence model of Baldwin and Lomax [102] are used, and the transition is fixed 
at 3% of chord. The overall agreement with experiments is quite good. The i/ 2 -norm residual of 10” 7 
can be reached in about 900 steps. 


5.6. A 3-D Steady-State Computations by a Point-Relaxation Implicit Method [103] 

Figure (5.25) shows pressure contours and energy contours in the plane of symmetry of the Aeroassist 
Flight Experiment (AFE) configuration at = 1429 m/s, p <*, = 60.136 N/m 2 , Too = 52.22°K, 
T w — 300°K, and Moo = 9.86. The configuration is a raked, elliptic cone with a circular shoulder. 
The body has a circular cross section when viewed perpendicular to the raked plane. The vehicle 
is designed as a test platform for a comprehensive series of experiments to define the flow field of 
an Aeroassisted Orbital Transfer Vehicle (AOTV) at high altitudes (above 75 km) returning from a 
Geosynchronous Earth Orbit (GEO). 

The numerical method employs a finite- volume, point-relaxation implicit procedure of the symmetric 
TVD formulation (limiter (4.33d)) of the governing Navier-Stokes equations. Gauss-Seidel iteration 
is employed across data planes in the sweep direction (from the body, accross the captured shock to 
the inflow boundary and back). Jacobi relaxation is used with respect to discretization within a single 
data plane. At each cell interface, Roe’s averaging is used to define eigenvectors and eigenvalues, and 
is set to a constant value of 0.2. Courant numbers up to 40 can be used to accelerate convergence. 
But Courant numbers of 1 to 2 must finally be used to damp high-frequency errors. The grid size 
(at close to the converged solution) was 64 cells between the body and the inflow boundary, 39 cells 
from the nose to shoulder, and 19 cells around the axis from 0° to 180°. The bow shock in this case 
is almost completely captured in two cells. 


5.7. The Entropy Condition and Blunt-Body Hypersonic Flows 

For Mach number ranges from 1.2 to 15, numerical experiments for one- and higher-dimensional 
unsteady flows containing unsteady shocks show that the second-order explicit TVD schemes via the 
local-characteristic approach are insensitive to the entropy correction for 0 < 6i < 0.1 in equation 
(3.18). In most cases 8i — 0 was used. For 0.1 < <!>! < 0.25, it has the effect of smearing the 
discontinuities slightly and sometimes has the advantage of improving stability in the sense of allowing 
a higher CFL number. 

For steady-state computations, 0.01 < 8i < 0.125 is a commonly used range for 0.6 < M < 4. The 
more complex formula suggested in reference [32] was also employed but shows no visible difference 
in shock resolutions or improvement in convergence rate when compared with the constant 8\ . 

For blunt-body steady-state flows with M > 4, the initial flow conditions at the wall are obtained 
using the known wall temperature in conjunction with pressures computed from a modified Newtonian 
expression. Also, for implicit methods, a slow startup procedure from free stream boundary conditions 
is necessary. Most importantly, it is advisable to use ^ as a function of the velocity and sound speed. 
In particular 


(My+i “ + K+|l + c j+|) ( 5 * 9a ) 

(^i)k+i “ ^(l u jfc+i| + l v *+i| + c fc+i) (5.9b) 

with 0.05 < 8 < 0.25 appears to be sufficient for the blunt-body flows for 4 < M < 25. Equation 
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(5.9) is written in Cartesian coordinates. In the case of generalized coordinates, the u and v should be 
replaced by the contravariant velocities, and one half of the sound speed would be from the ^-direction 
and the other half would be from the redirection. For implicit methods, it is very important to use 

(5.9) in ip(z) on both the implicit and explicit operators. For the implicit operator, instead of using 

(5.9) for t/>(z), numerical experiments showed that replacing the max/ t />( a y+i) and max/ t f>(a l k +±) 

in equations (5.8g) and (5.8h) by 


l Ht i \ + \ v i+h\ + c Hi 

(5.10a) 

k+{ \ + K+^l + C k+i 

(5.10b) 


performed well for high Mach number blunt-body flows. Again, in the case of generalized coordinates, 
the u,v and c in (5.10) should be the corresponding transformed coordinate variables. 

Numerical experiments (perfect and equilibrium real gases) [104,105] for either explicit or implicit 
methods showed that increasing the coefficient 8 has the effect of speeding up the convergence rate 
but smears the discontinuities slightly, since the bigger the 8 , the larger the numerical dissipations 
being added. One way of improving the convergence rate (without employing additional relaxation 
methods or multigrid procedures but retaining sharpness of the discontinuities) is to first use a larger 
8 in the range of 0.15 < 8 < 0.4 to obtain a converged solution and then use 0.01 < 8 < 0.15 as a 
postprocessor to obtain better resolution at discontinuities. Figures (5.26) and (5.27) show examples 
of two-dimensional steady-state blunt-body results. Figure (5.26) is a viscous computation and figure 
(5.27) is an inviscid, equilibrium, real-gas computation using (5.9). Details of the computations and 
related numerical examples can be found in references [104,105]. The following are short descriptions 
of the numerical results for figures (5.26) and (5.27). 

Figure (5.26) shows the Mach contours computed by the explicit MacCormack method and the 
implicit symmetric TVD scheme (5.8) using limiter (4.33c) of a viscous shock on shock interaction 
on a blunt cowl lip in the low hypersonic range. Extensive study on flow fields of this type were 
reported in references [106-108]. The computed result by shock-fitting for the outer bow shock and 
MacCormack methods for the interior of reference [107] is illustrated here for comparison. This flow 
field is typical of what will be experienced by the inlet cowl of the National Aerospace Plane (NASP). 
The free stream conditions for this flow field are M <*> = 4.6, Re^oo — 10,000, Pr = 0.72, p 0 0 = 14.93 
N/m, = 167°K, T w = 556°K, and 7 = 1.4 for a perfect gas. An oblique shock with an angle 
of 20.9° relative to the free stream impinges on the bow shock. Various types of interactions occur 
depending on where the impingment point is located on the bow shock. As shown by the Mach 
contours, the impinging shock has caused the stagnation point to be moved away from its undisturbed 
location at the symmetry line. The surface pressures and heat transfer rates at the new stagnation 
point can be several times larger than those at the undisturbed location of the stagnation point. In 
addition, a shear layer emanates from the bow shock and impinging shock intersection point and is 
intercepted by a shock wave which starts at the upper kink of the bow shock. The interacting shock 
waves and shear layers are confined to a very small region and must be captured accurately by the 
numerical scheme if the proper surface pressures and heat transfer rates are to be predicted correctly. 
The convergence rate using scheme (5.8) in this case is many times slower than the airfoil calculations 
discussed in section 5.4. The grid size for the MacCormack method is 51 X 31 and for the symmetric 
TVD is 151 X 71. The computational domain and number of grid points are larger for the symmetric 
TVD method than those for the MacCormack method, since the bow shock is also captured by the 
scheme. 
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Figure (5.27) shows the density contours computed by a second-order, non-MUSCL, upwind TVD 
scheme (4.40) using limiter (4.34c) of an inviscid blunt-body flow in the hypersonic equilibrium real 
gas range. The free stream conditions for this flow field are M 0 0 = 25, — 1.22 X 10 3 N/m 2 , 

Poo = 1.88 -2 kg/m 3 , and Too — 226°K. The grid size is 61 x 33 for the full (half) cylinder. The shock 
is at approximately fourteen points from the wall on the symmetry axis. The relaxation procedure 
employs a second-order Runge-Kutta method with a CFL of 0.5. The generalized Roe’s average [57] is 
used to define eigenvectors and eigenvalues, and 8 is set to a constant value of 0.15. Figure (5.27) shows 
the solution at 5000 steps with a four order of magnitude drop in the I^-norm residual. Pressure and 
Mach number contours converge and stabilize after 3000-4000 steps but the convergence rate is much 
slower for the density. The bow shock is captured in two to three grid points. Since the temperature 
of this problem can go up to 8100°K and the curve fits of Srinivasan et al. [89] used are accurate only 
for temperatures up to 6000° K, there is an uncertainty in the accuracy of the computed result and 
further investigation is needed. 

In summary, steady-state computations for higher Mach number, blunt-body flows in the hypersonic 
regime are very sensitive to initial conditions and numerical boundary treatments. A proper choice 
of the entropy parameter 8\ for higher Mach number flows not only helps in preventing nonphysical 
solutions but can act as a control of the convergence rate and helps in the resolution of discontinuities. 
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VI. Efficient Solution Procedures for Large Systems with Stiff Source Terms [99] 


In the application of modern shock-capturing methods like the TVD type of schemes to the chemi- 
cally reacting flows, Carofano [81] was the first to introduce the formalism that enabled full coupling 
in Marten’s explicit TVD scheme for a two-species, two-dimensional unsteady flow in Cartesian co- 
ordinates. However, due to the system size and the varying time-scale nature of the problem, the 
operations count increases nonlinearly as the number of species increases. To avoid solving a large 
system, Gnoffo and McCandless [109] and Gnoffo et al. [103] uncoupled the species equations from 
the fluid dynamics equations and solved these two sets of systems of nonlinear partial differential 
equations in a time-lag fashion (loosely coupled method) by using a point-relaxation technique with 
a second-order symmetric TVD scheme of Yee [10,87] and an upwind TVD scheme of Osher and 
Chakravarthy [11]. Eberhardt and Brown [110] attempted to use the eigenvalues and eigenvectors of 
the fluid dynamics equations alone to obtain a “fully coupled” first-order explicit TVD scheme for a 
one-dimensional flow. The results of Eberhardt and Brown showed excessive smearing at the shock 
when compared with the true, fully coupled explict TVD result. Their motivation for designing such a 
coupling procedure was to optimize the operations count by avoiding multiplication of large matrices. 
However, as will be demonstrated below, if one makes use of the unique structure of eigenvectors and 
eigenvalues for fluid flow of this type, the fully coupled formulation can be simplified even for a large 
number of species, thus providing a more efficient solution procedure than one might have anticipated. 
Moreover, using the eigenvalues and eigenvectors for the fully coupled equation set allows one to have 
the freedom of controlling the proper amount of numerical dissipation for the individual waves [10]. 
In particular, for the two-dimensional chemically reacting flows, the number of linear waves is ns + 1 
in each spatial direction where ns is the number of linearly independent species. Note that in or- 
der to capture contact discontinuities accurately, it is very important to apply the proper amount of 
numerical dissipation to the linear waves. 


Two types of schemes are proposed. If the stiffness is entirely dominated by the source term, a semi- 
implicit TVD type of shock-capturing method is proposed for steady-state calculations provided that 
the Jacobian of the source terms possesses certain properties. The proposed semi-implicit scheme can 
be viewed as a variant of the Bussing and Murman point-implicit predictor-corrector scheme [50] with 
a more appropriate numerical dissipation for the computation of strong shock waves in the hypersonic 
regime and a speed up in the convergence rate for steady-state applications. The predictor-corrector 
scheme of Bussing and Murman in turn is the explict MacCormack scheme with the source term 
treated implicitly. 


However, if the stiffness is not solely dominated by the source terms (e.g., stiffness due to highly 
irregular grid and/or viscous flows), a fully implicit method would be a more efficient procedure. The 
situation is complicated by problems in more than one space dimension, and the presence of stiff source 
terms further complicates the solution procedures for alternating direction implicit (ADI) methods. In 
fact, there seems to be no straightforward way of efficiently treating general stiff source terms implicitly 
with ADI procedures. Several alternatives will be discussed. The proposed fully implicit relaxation 
algorithm can be viewed as a variant of a fully coupled form of the algorithm proposed by Gnoffo and 
McCandless [109]. An implicit algorithm with explicit coupling between fluid and species equations 
proposed by the author will also be stressed here. Many existing perfect-gas or equilibrium real-gas 
computer codes can easily be modified to include this algorithm, which is a compromise between the 
loosely coupled implicit method of [109] and the fully coupled, fully implicit TVD method proposed 
here. To make this section more self-contained, some of the variables that were defined earlier will be 
repeated. 
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6.1. An Explicit Predictor-Corrector Algorithm for Systems with Source Terms 

In this section, all the necessary terms that are required for the basic TVD scheme for the com- 
pressible inviscid chemically reacting flow equations are derived. 

The Governing Equations : Consider a two-dimensional system of nonhomogeneous hyperbolic conser- 
vation laws, 


dU dF{U) dG(U) 
dt + dx ^ dy 


= S(U). 


( 6 . 1 ) 


Here U, F(U ), G(U), and S(U ) are column vectors of k components. Let A = dF/dU and B = 
dG/dU , with (a*, a 2 , ...,aj:) and (a y , a 2 , ..., a y ) being the eigenvalues of A and B. Denote R x and R y 
as the matrices whose columns are eigenvectors of A and B , and denote R~ x and R~ x as the inverses 
of R x and R y . In the case of the compressible inviscid flow equations with chemical reactions, the 
global continuity equation is replaced by the individual species continuity equations, 
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Here m = pu , n = pv , and s l represents the production of species from chemical reactions. The 
variables are the velocity components u and v, the pressure p, the total energy per unit volume e, and 
the density of the *th species p l . Also, p = an< ^ C V — P z ■> where ns is the number of species in 

the model and c x is the species mass fraction. Equation (6.2) assumes the p x are linearly independent. 


The eigenvalues of A and B are 


(<4,...,a”* +3 ) = (u,...,«,u + a,u,u - a), 
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Here the so called “frozen sound speed” a is 

a 2 = p p Ap e {H - u 2 - v 2 ), 


with 


E ns . 

i=l CPp ' 

E fl s 

c l = 1, 
dp 

** = V 
dpi 




Pe = 


de 


(6.3a) 
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m,n,p 



(6.9) 


e + p 
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The frozen sound speed a defined in (6.4) has no physical meaning. It is defined here for the convenience 
of notation for the basic scheme. 

The superscript n is used for the time-index and should not be confused with the n = pu in equation 
.2). Let a 1 t , iL , i , R~ l v denote the quantities a 1 R x , R~ l evaluated at some symmetric average 

> 3+2 2 3+2 

of U Jik and Uj+ i,fc. Similarly, let a l k+l , R k+ i, denote the quantities a l yi R y , R~ x evaluated at 

2 2 • _ . 
some symmetric average of Uj,k and U ■ In the case of chemically reacting flows, Rj+i, 

R- 1 , a 1 , , Rl. , i , and R7 1 1 are defined in ways similar to the ones used by Huang [30] and Carofano 

3+ i k + 2 2 k +2 

[81]. Also, although the commonly used Roe’s average is no longer valid for a nonperfect gas, a close 
approximation based on a generalized Roe’s average for equilibrium flows as discussed in section 4.3 
and reference [57] can also be used. 

For thermally and chemically nonequilibrium flows, the eigenvalues and eigenvectors have a similar 
structure. For the two-dimensional system (6.1), if nt is the number of thermal energy variables, then 
the eigenvalues in the z-direction will have (ns + 1) “u” characteristics plus u + a and u - a 

characteristics. Here the values “a” will reflect the added thermal energy variables 

As discussed in sections 3.3, 3.6, and 3.7, the TVD property is only valid for homogeneous scalar 
hyperbolic conservation laws. Certain types of source terms might preserve the original TVD property 
of the homogeneous PDE, and others might not. However, disregarding the type of bounded source 
terms, one is not precluded from the application of TVD schemes when source terms are present. 
But, extreme precaution has to be taken (especially for explicit TVD schemes) in the procedure of 
including the source terms. This applies particularly to stiff source terms of the types in thermally 
and chemically nonequilibrium flows. Procedures for time- accurate computations for problems with 
stiff source terms are currently under investigation. However, in the use of time-accurate methods as 
relaxation procedures for steady-state computations, the semi-implicit algorithm discussed in section 
3.7 appears to work quite well. 

Explicit Preditor- Corrector TVD Scheme : With the above notation, a formal extension of the scalar 
explicit second-order TVD method (section 3.6) in predictor-corrector form via the local-characteristic 
approach for the nonlinear hyperbolic system (6.1) with nonzero source terms can be written as 
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(6.10b) 

(6.10c) 

(6.10d) 

(6.10e) 


The elements of the vector 3> J+ i are the same as equations (3.59d) and (3.59e) with a J+ ± replaced 
by a l 3+ , Aj + x replaced by a*. + i , and Qj+x replaced by Q l .+ k . Here Ay is the grid spacing such 
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that yk = kAy. Research is underway to study the influence of the stiff source terms on the TVD 
correction step for the second and third terms of (6.10e) to be evaluated at u n compared with (6.10e). 
For steady-state calculations and for a contractive type of source term S(U), it appears that (6.10e) 
might helped improve the convergence rate by using the most updated information in the relaxation 
procedure. Other procedures are currently under investigation. 

6.2. More Efficient Solution Procedures for Large Systems 

The extra computation in (6.10) compared with a classical central-difference shock-capturing scheme 
such as the Lax-Wendroff method is due to the vectors (R^) J± a. At first glance, the vectors 
and i involve matrix and vector multiplication of dimension ns + 3 for every grid point, and 

thus tend to discourage their adoption in problems other than ideal gas flows. Researchers such as 
Gnoffo and McCandless [109], and Eberhardt and Brown [110] were motivated to pursue other avenues 
to solve the complicated chemically reacting flow problems. However, as will be demonstrated in this 
subsection, if one makes use of the unique structure of the eigenvectors and eigenvalues for fluid flow 
of this type, the fully coupled formulation can be simplified even for a large number of species, and 
thus becomes a viable approach. 

With straightforward manipulations, the computation for scheme (6.10) can be simplified tremen- 
dously. The corresponding vector a in equation (4.12c) for system (6.1)-(6.2), for example, can be 
expressed as 
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Here, for example, (Ap l ) J+ ± k = - p) )k , and it is understood that the spatial indices in (6.11)- 

(6.13) are at (j + 2 , k). Similarity, also has a very simple form. For instance, the R& associated 

with the x-direction flux can be expressed as 
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Here the spatial indices on (6.14) are at (j + As one can 

see, the terms in equations (6.11) and 


(6.14) due to the species equations are simple and do not require many operations. Therefore, the 
increase in the number of species equations is not as CPU-intensive as one might have anticipated. 


6.3. A Semi-implicit Predictor-Corrector TVD Algorithm and a 3-D Example 


As mentioned in discussing the scalar schemes, the explicit TVD scheme (6.10) can be used for either 
time-accurate or steady-state calculation. It is second-order accurate in time and space. However, if 
the source term is stiff, the restriction in the time-step due to stability requirements is prohibitively 
small and (6.10) is not practical, especially for steady-state applications where a large system such 
as (6.2) is involved. In this section, a semi-implicit method is proposed. Another alternative is a 
fully implicit method. The basic implicit scheme and the related difficulty in efficiently extending the 
implicit method to two dimensions with stiff source terms will be discussed in the next section. 

If one follows the idea of Bussing and Murman [50] in treating the source term implicitly, a semi- 
implicit predictor-corrector TVD algorithm for steady-state computations can easily be obtained. It 
can be written as a one-parameter family of time-differencing schemes for the source term; i.e., the 
following formulation includes scheme (6.10). The proposed scheme for system (6.1) can be written as 
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Here, D is assumed to be invertible; i.e., only source terms with Jacobians such that D is invertible at 
each grid point are permissible. Again, 9 has the same meaning as in the scalar case. For 9 ^ 0, the 
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source terms are treated implicitly. If 0 — 1, the time-differencing for the source term is first-order and 
(6.15) is best suited for steady-state calculations. See sections 3.3, 3.6, 3.7 and 6.1 for the discussion 
of equation (6.15d) and the usage of the scheme (6.15) in general. 

One can simplify equation (6.15) by partitioning the vectors U,F,G,S , and D in equation (6.2) as 
follows: 
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Here D 21 is a null matrix and D 22 is an identity matrix. With the above definitions, the scheme can 
be greatly simplified. The procedures are as follows: taking the predictor step, for example one first 
solves for (At/ 77 )* 1 ) by 
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where ( r.h.s) 7 is the right-hand side of (6.17a) with all the indices “H” replaced by “I”, and with the 
term At(S* k ) n added. In other words, one only has to invert the D 11 matrix of dimension (ns, ns) 
instead of (ns + 3, ns + 3). Similarly, one can simplify the corrector step in the same way. The solution 
obtained from the above procedure is used in (6.15d). Or, to explain it in another way, one solves the 
predictor step of the fluid equations 


dU u dF n (U) dG n (U) _ 

dt dx ^ dy 

explicitly, then uses the result to solve the predictor step of the species equations 


(6.18a) 


d Jjy_ dF J (U ) dG’(U) _ j 
dt dx dy 


(6.18b) 


explicitly, with the exception that the chemical reaction terms are treated implicitly in (6.18b). One 
then repeats the same procedure for the first corrector step. The solution obtained from the first 
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corrector step is then used to solve for a and #$ in equations (6.11)-(6.14) for the complete system 
(6.1) so that one can solve for the second corrector step. Here, the first corrector step means the step 
to obtain from U^ l \ and the second corrector step means the step to obtain U n+1 from in 
(6.15). Note that the second corrector step is the important part of the algorithm that deviates from 
the Bussing and Murman method. This step, which ensures that the method will have the TVD-type 
properties, is designed to capture shock waves without the associated spurious oscillations. 

In the case where S 11 is not a null vector and it is not stiff, equation (6.17) is still applicable, except 
one has to add the term At(S^ k ) n to the right-hand side of (6.17a). In steady-state calculations 
where body-fitted coordinates are used, one can further speed up the convergence rate by using a local 
time-stepping approach [18,101]. 

To verify the current approach and to make a fair comparison with a known method, the proposed 
semi-implicit scheme was implemented into an existing three-dimensional code [111]. This existing full 
Navier-Stokes code, originally developed by A. Kumar at NASA Langley Research Center, contains 
the explicit MacCormack scheme with source terms treated implicitly. A detailed description of a 
numerical experiment and the extension of (6.15) to three-dimensional, chemically reacting flows in 
generalized coordinates is given in Shinn et al. [100]. Only partial results from reference [100] are 
presented here to illustrate the performance of the proposed semi-implicit scheme. In reference [100], 
three symmetric averages for the eigenvalues and eigenvectors were studied. Aside from a slight 
difference in convergence history, no visible differences in resolution were observed among the different 
averages. The computations shown employed =0.1 for the entropy correction function ip(z) in 
equation (3.18). 

Figures (6.1)-(6.3) show a preliminary test result for a three-dimensional, five-species viscous re- 
acting flow using the semi-implicit TVD scheme (6.15) together with (6.13f), as compared with an 
existing classical shock-capturing method which supplies numerical dissipation linearly [111]. The 
numerical result for the semi-implicit TVD method is shown to be oscillation-free around the shock, 
while the time spent per iteration is approximately double when compared with the method used in 
[111]. The configuration of the numerical experiment is shown in figure (6.1). Although this is a 
two-dimensional flow, to check out the three-dimensional code it is computed as a three-dimensional 
flow with the appropriate boundary conditions in the y-direction. A uniform grid consisting of 31 
points in the ^-direction, 6 points in the y-direction and a viscous grid consisting of 51 points in the 
^-direction was used. The inflow conditions are the pressure p = 1 atm, the temperature T = 1200 K, 
and the Mach number M — 4 for a premixed air and hydrogen fuel. The species considered are H 2 , 
02, OH , # 2 0, and iV 2 , with the two reactions (iV 2 being inert) 


#2 + 0 2 ^ 2 OH, 
20 # + # 2 — 2 # 2 0 . 
The reaction rates for the above are given in reference [111]. 


6.4. A Fully Implicit TVD Method and a 3-D example 

Another type of shock-capturing scheme that might be applicable to (6.1) is a one-parameter family 
of explicit and implicit TVD-type of schemes (section 3). For the nonequilibrium equation (6.1), these 
schemes can be written as 
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(6.19) 
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Here 9 has the same meaning as before. The numerical fluxes F J± i k and G jk± i have the same 
meaning and form as the scheme for the homogeneous PDE’s (5.1). Again, the MUSCL or non- 
MUSCL approaches as discussed in section 4.5 for the implicit method are applicable here. 


6.4.1 A Conservative Linearized Form For Steady-State Applications 

To solve for U n+1 in (6.19), one normally needs to solve a set of nonlinear algebraic equations 
iteratively. One way to avoid this is to linearize the implicit operator, and solve the linearized form 
by other means. Following the same procedure as in section (3.8) and in Yee [10,43], a conservative 
linearized form of (6.19) can be written as 
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(6.20a) 


E = U n+1 - U n , (6.20b) 

where the F j+±,k anc * G j,k+± have the same form as in equation (5.8) and section 

(4.5) except all the vectors are in Cartesian coordinates. 


6.4.2. Stiff Source Terms, ADI Approaches and Relaxation Methods 

The stiff terms D* k on the implicit operator (6.20) complicate the solution procedures for the 
commonly used ADI procedures. Normally, if D^ k are not stiff, one can reformulate (6.20) by an ADI 
procedure like the Beam and Warming [17] algorithm for an efficient solution process. Unfortunately, 
the Dj,k considered here is stiff; consequently the additional higher order terms due to the ADI 
formulation can no longer be ignored. In a different context, Van Dalsem and Steger [112] suggested a 
remedy if Dj k is a diagonal matrix with identical diagonal elements. However, for chemically reacting 
flows, the matrix D 3 k is full for the upper (ns, ns + 3) entities and no straightforward efficient way of 
utilizing ADI approaches for nonlinear system cases with a general stiff source terms can be found. 

The straightforward way of iteratively solving (6.19) as a set of nonlinear system of equations 
iteratively, or solving the linearized form (6.20) or similar non- ADI forms, appears to be quite expensive 
for large systems. Recently, Gnoffo et al. [103] successfully demonstrated the usefulness of a point- 
relaxation method on the implicit symmetric TVD scheme similar to (6.20) for a loosely coupled 
chemically nonequilibrium flow. Here a similar point-relaxation or line-relaxation method is proposed 
for the fully coupled system (6.15). Despite the fact that a larger equation set is involved than in 
[103], the extra operations are minimized by making use of the simplification procedure of section 
6.2. For a point-relaxation method, the size of the matrix inversion for (6.20) is (ns + 3, ns + 3) as 
opposed to the loosely coupled method of [103], where the size is (ns, ns). The gain in the freedom 
of controlling the appropriate amount of numerical dissipation for each individual wave more than 
compensates for the extra expense. More importantly, solving the fully coupled system is believed to 
have a better convergence rate than the loosely coupled approach. All of the necessary terms required 
for the implicit_scheme (6.20) are derived in section II. The implicit operator of (6.20a) is diagonally 
dominant for Dj^ = 0. Therefore, one has to make sure that the source term does not destroy 
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the diagonally dominant property which is required for some relaxation methods. In the following 
subsection, an even simpler form than (6.20) is proposed. 

6.4.3. An Implicit Algorithm with Explicit Coupling between Fluid and Species Equa- 
tions 

To avoid the inversion of large matrices, one can further simplify (6.20) by requiring that the coupling 
between the fluid and species equations be explicit. With this relaxed requirement, one effectively 
solves the fluid and species equations separately. Unlike the loosely coupled method used in [103] 
or the chemistry-split technique used in [113], the eigenvalues and eigenvectors of the fully-coupled 
equations are coupled explicitly between the fluid and species equations. This can be accomplished by 
partitioning U,F,G,S,D in the same way as in (6.16) and partitioning the Jacobians A and J3, and 
the numerical fluxes F j± l k and G - k± ± similarly. For example, the Jacobian A and the numerical 

fluxes F J± i k can be partitioned as 
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The dimensions of the subvectors and submatrices are the same as (6.16) if equation set (6.2) is chosen, 


and the A 11 , A 12 , A 21 , and A 22 

in this case are 
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Similarly, E 3> k and Ct x can be partitioned as 


(6.23c) 


(6.23d) 
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(6.24) 
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n x = 


(n x ) n (n x ) 


(rr ) 21 (n x ) 


12 


22 


(6.25) 


Recall that two choices of equation sets are available: namely, equation set (6.2), or the full fluid 
equations (i.e., keeping the global continuity equation) plus the ns - 1 species equations. The procedure 
suggested here is best suited for the latter equation set, since modification of many existing implicit 
solver computer codes (perfect gases or equilibrium real gases) would be minimal. The dimensions of 
subvectors and submatrices for the latter equation set would be slightly different. For example, the 
dimensions for A 11 , A 12 , A 21 , and A 22 would be (ns - 1) x (ns - 1), (ns - 1) x 4, 4 x (ns - 1), 

and 4 x 4, respectively. 

The procedures for either equation set are as follows. One solves the fluid equations (e.g., (6.18a)) 
implicitly, and then uses the result (U u ) n+1 to solve the species equation (e.g., (6.18b)) implicitly. In 
other words, one solves ( U n ) n+1 from the following: 





= (rhs) 11 7 (6.26a) 
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Ax 
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(6.26b) 
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The explicit coupling between the fluid and species equations is imbedded in (6.26b). After obtaining 
(U 11 )^ from (6.26), the solution is used to solve the species equation as follows: 


with 


I + 


~ - (*W 7 ) 

= (rhs) 1 - ( Ihs ) 12 , 


- M ^) 11 


E 1 

(6.27a) 
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+ fK^+i ) 12 " (n h + d n - + HVP 12 ) 


rn . 12 


~0At(D j>k ) \E U (6.27b) 


and 


(6.27c) 

The term (rhs) 1 is the right-hand side of (6.26b) with all the indices “II” replaced by “I” and the source 
term At(5^ fc ) 7 added. Here the quantities in equation (6.27) are evaluated with updated information 
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from equation (6.26). To simplify the procedure even further, one can drop the second term on the 
right-hand side of (6.27a) entirely; i.e., one can set equation (6.27b) to zero. In other words one can 
solve for ([/ J ) n+1 without using the most updated information from (6.26a); thus E 11 in (6.27b) in 
this case is zero. 

With this simplified procedure, one only has to solve two reduced systems of dimension 3 and ns 
(or 4 and ns - 1). For line-relaxation methods and ns > 3, this procedure can provide a large savings 
in operations count. Although numerical experiments on this simplified procedure have not yet been 
done, one would expect that the current method will give a faster convergence rate than the method 
of Gnoffo et al. [103] which was demonstrated to be applicable to many three-dimensional blunt- 
body problems. A numerical computation by Gnoffo et al. [103] will be presented in section 6.4.4. If 
point relaxation were used, the only difference between the two methods is that reference [103] uses the 
eigenvalues and eigenvectors of the individual subsystems (fluid and species) alone, whereas the current 
method uses the eigenvalues and eigenvectors of the full system. The use of the full eigenvalues and 
eigenvectors set for the current method is believed to enhance the coupling between the two systems 
without imposing additional conditions as employed in reference [103]. 

For steady-state application, an algorithm utilizing the TVD scheme for viscous flows is to difference 
the hyperbolic terms the same way as before, and then central difference the viscous terms. The final 
algorithm is the same as equation (6.20) (or (6.26)-(6.25)),except that the spatial central differencing of 
the viscous term is added to the right-hand side of equation (6.20) (or (6.26)-(6.27)). Numerical tests, 
comparison with other approaches, and recommendations will be reported in a future publication. 


6.4.4. A Numerical Example for a Loosely Coupled Point-Relaxation Implicit Method [103] 

To illustrate the applicability of a point-relaxation implicit algorithm using the symmetric TVD 
scheme (4.33) for three-dimensional chemically nonequilibrium flows, Gnoffo et al. applied the scheme 
to solve a configuration similar to the ones used in section 5.6. However, a free stream Mach number 
of 32 was used in this computation. The free stream conditions for this case were = 8917 m/sec, 
Poo = 1.54 N/m 2 , and T = 197°X. The maximum body diameter was 14 ft., corresponding to the 
full-scale Aeroassist Flight Experiment (AFE) vehicle. The kinetic model of Dunn and Kang [114] was 
used which involves 11 species, (N,0, ^2,02, NO, N + ,0 + , iV^,0^, NO + , e) and 26 reactions. 

Figure (6.4) shows contour plots of electron number density. The global contour plot in figure (6.4a) 
serves only to define the shock layer. The blowup of the shock layer near the far shoulder in figure 
(6.4b) shows that contour lines run nearly parallel with the body and shock. Number densities vary 
from 2 x 10 13 /cm 3 at the body to 10 15 /cm 3 at the shock. Initial conditions for this nonequilibrium 
chemistry test case were taken from a converged fine-grid perfect-gas solution. 

The same relaxation procedure and limiter as described in section (5.6) were used, except now 
the species equations are solved separately from the fluid equations. Roe’s averaging is used to 
define eigenvectors and eigenvalues for the fluid and species subset of equations separately (i.e. the 
full eigenvectors and eigenvalues of the fully coupled system were not used). Other thermodynamic 
derivatives’ average values were defined in reference [103]. Procedures which allow a good control of 
this loosely coupled procedure were discussed in reference [103], 
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VII. Concluding Remarks 


A unified formulation accompanied by practical fluid dynamics computations of a class of shock- 
capturing methods for one- and higher-dimensional nonlinear systems has been presented. In time- 
accurate computations where experimental data were available, good agreement between the numerical 
results and the experimental data were observed. Steady-state computations performed well for sub- 
sonic, transonic and supersonic flows. However, the convergence rate for hypersonic flows is less than 
desired and improvement in this area is the current pacing item. 

The applicability of these shock-capturing methods for equilibrium real gases was studied. Prelim- 
inary results for one-dimensional shock-tube and two-dimensional steady-state blunt body problems 
show that the shock and contact discontinuity resolution were not affected by the state equation for 
a wide range of flow conditions. 

Two numerical algorithms for hyperbolic conservation laws that are suitable for large systems of 
thermally and chemically nonequilibrium steady-state flows in the hypersonic regime were proposed. 
The specific properties of the governing equations for fluid flow of this type were taken into consid- 
eration for more efficient solution procedures. The main areas of consideration were to minimize the 
operations count, increase the allowable time-step constraint imposed by the stiff source terms, and 
expand the shock-capturing capability beyond classical approaches. Details of all these considerations 
were described. A preliminary test problem showed certain advantages of the proposed semi-implicit, 
high-resolution shock-capturing scheme over the classical ways of supplying numerical dissipation. 
More numerical testing and study will be pursued in the immediate future. 

In summary, the performance of the schemes presented for one- and two-dimensional gas-dynamics 
problems in conjunction with the various Riemann solvers can be divided into the following aspects. 

Looking at the Riemann solvers, in general the advantages of the generalized van Leer flux-vector 
splitting over the generalized Steger and Warming formulation remain for a real gas, with slightly less 
dissipation at the discontinuities. The local-characteristic approach (approximate Riemann solver) 
gives results very similar to the generalized van Leer flux-vector splitting formulation. 

For one-dimensional problems, the difference in computational effort required by the three Riemann 
solvers is small. The main difference lies in the MUSCL and non-MUSCL approaches. The operations 
count between the non-MUSCL and MUSCL is within 30% for a perfect gas. However, due to extra 
evaluation in the curve fitting between the left and right states in an equilibrium real gas for the 
MUSCL formulation, additional computation is required for the MUSCL approach. The amount of 
additional computation increases nonlinearly as the spatial dimension increases. 

Looking at the numerical schemes, the main difference seems to occur between the upwind and the 
symmetric approaches. The upwind schemes give better results for contact discontinuities. On the 
other hand, symmetric schemes have a better stability and produce shock resoltuion similar to that of 
the upwind schemes; yet they require less CPU time per time-step and are less sensitive to numerical 
boundary condition treatments. 

In the case of a non-MUSCL approach, limiter (4.33e) is the most accurate among (4.33c) - (4.33e) for 
the symmetric TVD scheme. As for the upwind schemes, limiters (4.34d) and (4.34e) are very similar, 
whereas limiter (4.34f) gives very accurate results for contact discontinuities but is too compressive, 
causing slight oscillations in smooth regions for high Mach number cases. A combination of limiters 
such as (4.34d) or 4.34e) for the nonlinear fields and (4.34f) for the linear field seems to be a good 
compromise. In the case of the MUSCL approach, only limiters (4.34c) and (4.34e) were studied. 
Between the two limiters, (4,34e) produces higher shock-resolution than (4.34c). 
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None of the differences for the various approaches observed for the explicit schemes seems to be 
decisive for the one-dimensional tests, but factors such as stability or computational efficiency need 
further investigation in multidimensional tests. The main differences between the methods lie in 
their versatility in extending to implicit methods with efficient solution procedures, especially for 
multidimensional steady-state computations. Preliminary study shows certain advantages of the local- 
characteristic approach over the flux-vector splitting approaches. 

There is an important distinction between the flux-vector splittings and the local-characteristic 
approach for implicit methods. Unlike flux-vector splitting approaches, implicit methods employing 
the local-characteristic approach (non-MUSCL or MUSCL with first-order implicit operator such as 
(4.41f)) do not require the Jacobian of the fluxes. In many instances, the Jacobian of is 
relatively difficult to obtain. A similar difficulty applies to the MUSCL formulation via the exact 
Riemann solver or local-characteristic approach (if a second-order implicit operator is desired). 

Another important fact is that flux-vector splittings make use of the sound speed only, whereas 
local-characteristic approach of the Roe-type make use of the thermodynamic dervatives x and K °f 
equation (4.2). These thermodynamic derivatives put more stringent requirements on the curve fit 
that represents the thermodynamic properties of the gas. In this regard, the curve fits of Srinivasan 
et al. may be deficient for the approximate Riemann solver as can be seen from figure (4.1), case D. 
One probably needs more improved curve fits than those of reference [89] before a definite conclusion 
can be drawn about the accuracy of the different Riemann solvers and schemes for real gases. 
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Fig. 3.3 Graphical representation of three limiters (3 52,b,c,d). 
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Fig. 4.2 Comparison of the schemes for test case B, for a perfect gas. 

(a) Symmetric TVD, non-MUSCL; (b) upwind TVD non-Ml SC'L. (c) upwind TVD, MUSCL 
(d) van Leer splitting, MUSCL; (e) Steger and Warming splitting. MUSCL. 

Note: (a)-(c) use the identical approximate Riemann solver. 




Fig. 4.3 Comparison of the schemes for test case E, for a perfect gas. 

(a) Symmetric TVD, non-MUSCL; (b) upwind TVD, non-MUSCL; (c) upwind TVD, MUSCL; 
(d) van Leer splitting, MUSCL; (e) Steger and Warming splitting, MUSCL. 

Note: (a)-(c) use the identical approximate Riemann solver. 
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Fig. 4.4 Comparison of the schemes for test case C, for a perfect gas. 

(a) Symmetric TVD, non-MUSCL; (b) upwind TVD, non-MUSCL; (c) upwind TVD, MUSCL; 
(d) van Leer splitting, MUSCL; (e) Steger and Warming splitting, MUSCL. 

Note: (a)-(c) use the identical approximate Riemann solver. 
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Fig. 4.5 Comparison of the schemes for test case B, for equilibrium air. 

(a) Symmetric TVD, non-MUSCL; (b) upwind TVD, non-MUSCL; (c) upwind TVD, MUSCL; 
(d) generalized van Leer splitting, MUSCL; (e) generalized Steger and Warming splitting, MUSCL 
Note: (a)-(c) use the identical approximate Riemann solver. 




Fig. 4.6 Comparison of the schemes for test case E, for equilibrium air. 

(a) Symmetric TVD, non-MUSCL; (b) upwind TVD, non-MUSCL; (c) upwind TVD, MUSCL; 
(d) generalized van Leer splitting, MUSCL, (e) generalized Steger and Warming splitting, MUSCL. 
Note: (a)-(c) use the identical approximate Riemann solver. 
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Fig. 4.7 Comparison of the schemes for test case C, for equilibrium air. 

(a) Symmetric TVD, non-MUSCL; (b) upwind TVD, non-MUSCL; (c) upwind TVD, MUSCL; 
(d) generalized van Leer splitting, MUSCL; (e) generalized Steger and Warming splitting, MUSCL 
Note: (a)-(c) use the identical approximate Riemann solver. 
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Fig. 4.8 Comparison of the schemes and limiters for test case B for a perfect gas. 
(a) First-order Roe; (b) second-order Harten (3.45), non-MUSCL; 

(c)-(g) upwind TVD (4.34), non-MUSCL, limiters (4.34c,d,e,f); 

(h)-(j) symmetric TVD (4.33), non-MUSCL, limiters (4.33c,d,e); 

(k) upwind TVD, MUSCL; (1) van Leer splitting, MUSCL; 

(m) Steger and Warming splitting, MUSCL. 

Note: (a)-(k) us Roe’s average, and (k)-(m) use the same limiter (4.34e). 
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Fig. 4.9 Comparison of the schemes and limiters for test case B for a real gas. x 

(a) First-order Roe; (b) second-order Harten (3.45), non-MUSCL; 

(c)-(g) upwind TVD (4.34), non-MUSCL, limiters (4.34c,d,e,f); 

(h)-(j) symmetric TVD (4.33), non-MUSCL, limiters (4.33c,d,e); 

(k) upwind TVD, MUSCL; (1) generalized van Leer splitting, MUSCL; 

(m) generalized Steger and Warming splitting, MUSCL. 

Note: (a)-(k) use the identical approximate Riemann solver, and (k)-(m) use the same limiter (4.34e). 
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(a) Explicit TVD method. 


(b) Implicit TVD method. 



(c) Conventional implicit method 


(d) First-order implicit flux- vector 
splitting method. 


Fig. 4.10 Density distribution of a quasi-one-dimensional divergent nozzle for four different methods. 
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Fig. 5.1 Schematic of a shock wave diffraction problem. 
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Fig. 5.2 The 317 x 100 orthogonal grid. 
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Fig. 5.3 Density contours computed by a time-split, non-MUSCL second-order 
upwind TVD scheme. 
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Fig. 5.4 Density contours computed by a predictor-corrector symmetric TVD scheme. 
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Fig. 5.5 Comparison of density contour plots from different TVD schemes. 

(a) Predictor-corrector symmetric TVD with limiter (4.33c); 

(b) predictor-corrector symmetric TVD with limiter (4.33d); 

(c) predictor-corrector symmetric TVD with limiter (4.33e); 

(d) time-split second-order upwind TVD, non-MUSCL with limiter (4.34d,f) 
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Fig. 5.7 Density contours (a) and pressure contours (c) computed by a predictor-corrector 
symmetric TVD scheme compared with the Schlieren photograph (b). 

Notes: T.P. — triple point; M.S. — Mach stem, V — vortex, C.D. — contact discontinuity 
Ft S. — reflected shock, I.C. = incident shock. 
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Fig. 5.8 Schematic of the experiment. 
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Fig. 5.9 Schematic of the shock tube computational domain (Cartesian grid). 
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Fig. 5.10 Schematic of the computational domain (curved incident shock). 
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Fig. 5.12 Density contours computed by the time-split upwind TVD scheme (4.34a,d,f) at six 
sequential stages of the diffraction process. 
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Fig. 5.13 Interferogram (a), Density contours (b), pressure contours (c), and grid (d) at an 
instant where a vortex is being formed near the trailing edge of the airfoil. 
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Fig. 5.19 Pressure contours (lei 
symmetric TVD sche] 
(a) Limiter (4.33c); (1 
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Fig. 5.20 Comparison of a symmetric TVD (SYMTVD) scheme with an upwind TVD (UPTVD) scheme for the 
NACA0012 airfoil with — 0.8, a — 1.25 using a 163 x 49 C grid. 
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Fig. 5.21 Comparison of a symmetric TVD (SYMTVD) scheme with an upwind TVD (UPTVD) scheme for the 
RAE2822 airfoil with M 0 0 = 0.75, a = 3 using a 163 x 49 C grid. 
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Fig. 5.22 Comparison of a symmetric TVD (SYMTVD) scheme with ARC2D (version 150) for 
the Mach contours, pressure contours and entropy contours of the NACA 0012 airfoil 
with Moo = 1.2, a = 7 using a 163 x 49 C grid as shown in figure 5.20. 















Fig. 5.23 Comparison of a symmetric TVD (SYMTVD) scheme with ARC2D (version 150) for 
the Mach contours, pressure contours and entropy contours of the NACA 0012 airfoil 
with Mqo = 1.8,0; — 7 using a 163 x 49 C grid as shown in figure 5.20. 
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(a) Explicit MacCormack (bow-shock - shock-fitting). 


(b) Implicit 2nd-order symmetric TVD (bow-shock - shock-capturing). 


Fig. 5.26 Two-dimensional viscous blunt-body flows (Mach contours) for a 20.9° shock impingement 
(type III-IV) for a free-stream Mach number of 4.6 and Reynolds number of 10,000. 
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Fig. 6.2 Pressure profiles along z — 0.29 and 2 = 0.42 cm. 
— Semi-implicit TVD method 
- - Classical shock-capturing method 
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